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ABSTRACT  The  objective  of  this  research  was 
to  develop  a  composite  finite  element  analysis  that  is 
capable  of  including  the  nonlinear,  elastic-plastic  re¬ 
sponse  of  the  bond  between  reinforcing  fibers  or  bars 
and  matrix  material.  This  bond  is  very  important  in 
determining  the  behavior  of  reinforced  weak  or  brittle 
matrix  materials,  such  as  reinforced-concrete  or  ce¬ 
ramic-matrix  composites.  A  composite  element  cap¬ 
tures  the  response  of  a  reinforced  material  without 
the  expense  of  a  discrete  model  which  treats  each 
constituent  with  different  element  or  material  types. 
This  report  documents  the  development,  implemen¬ 
tation,  and  verification  of  the  composite  element  with 
bond  modeling  capabilities. 

The  composite  continuum  element  employs  a  bond 
model  based  on  classical  non-associative,  elasticity- 
plasticity  theory.  The  element  is  two-dimensional, 
plane  strain,  with  reinforcement  in  one  coordinate 
direction.  The  element  has  two  added  degrees  of 
freedom  per  node:  bond  slip  and  dilation.  These 
local  deformations  are  approximated  as  continuous 
field  variables  in  the  finite  element  scheme. 

The  finite  element  scheme  requires  special  com¬ 
posite  material  constitutive  relations  which  include 
bond  slip  and  dilation  effects.  A  unit  cell  micro- 
mechanics  analysis  based  on  virtual  work  concepts 
was  used  to  derive  the  required  relations. 


The  inelastic  bond  constitutive  relations  were  evalu¬ 
ated  using  a  Reduced-Newton  algorithm  developed 
especially  for  this  bond  model.  The  global  finite  el¬ 
ement  solution  uses  an  incremental  Newton-Raphson 
process  to  solve  the  resulting  nonlinear  system  of  equa¬ 
tions. 

The  element  was  tested  at  various  stages  of  its 
development.  Specializing  the  bond  model  for  the 
case  of  perfect  bond  conditions  allowed  for  simple 
verification  of  the  underlying  micro-mechanics  rep¬ 
resentation  of  the  composite  system.  The  numerical 
implementation  of  the  plasticity-based  bond  model 
during  testing  reproduced  the  salient  features  of  the 
theoretical  bond  model  as  a  stand-alone  unit  prior  to 
its  implementation  in  the  finite  element  setting.  Af¬ 
ter  the  composite  finite  element  was  incorporated  in 
a  finite  element  code,  analyses  of  published  beam 
and  bond  tension  tests  were  conducted.  The  pre¬ 
dicted  results  compared  favorably  with  the  published 
data,  considering  the  limitations  of  modeling  the  be¬ 
havior  of  the  matrix  and  the  reinforcement  as  linear 
elastic.  The  beam  analyses  illustrated  the  ease  with 
which  the  composite  analysis  can  be  applied  to  prac¬ 
tical  problems.  The  bond  tension  tests  clearly  dem¬ 
onstrated  the  element’s  capacity  to  model  bond  be¬ 
havior. 
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1.  INTRODUCTION 


This  project  as  originally  conceived  was  to  span  5  years;  however,  early  in  the  project, 
changes  in  funding  priorities  required  a  substantial  shortening  of  the  proposed  research.  The  chief 
goal  of  the  original  project  was  the  numerical  implementation,  in  a  composite  finite  element 
analyses,  of  a  plasticity  bond  model  for  reinforced  materials.  The  model  was  to  include  both  bond 
slippage  and  dilation  and  was  to  be  of  the  general  form  of  the  model  developed  by  Herrmarm  and 
Cox  (1992)  for  reinforced  concrete.  The  applications  of  the  analysis  were  to  include  both 
reinforced  concrete  and  fiber-reinforced  composites.  When  the  overall  duration  and  magnitiidft  of 
the  project  were  reduced,  it  was  necessary  to  substantially  reduce  the  goals  of  the  project.  As  a 
result  of  this  redefinition,  the  study  was  restricted  to  circumstances  where  the  cyclic  plasticity 
features  of  the  model  could  be  ignored.  In  addition,  the  matrix  and  reinforcing  materials  were 
assumed  to  be  linear  elastic,  and  the  analysis  was  restricted  to  plane  strain  conditions.  Finally, 
examples  were  limited  to  simple  reinforced-concrete  tests  and  structures. 

While  the  above  reduction  in  the  overall  goals  of  the  project  limits  its  immediate 
applicability,  it  does  not  really  affect  the  primary  objective  of  demonstrating  the  feasibility  of 
incorporating  a  plasticity-based  bond  law  that  includes  bond  dilation  into  a  composite  finite  element 
analysis  for  reinforced  structures.  The  possible  future  extension  of  the  analysis  to  three- 
dimensional  problems  and  its  incorporation  into  the  analysis  of  inelasticity  material  models  for 
matrix  and  reinforcmg  materials  are  all  steps  that  are  well  understood  (albeit  rather  long  and 
tedious).  Finally,  the  application  to  fiber-reinforced  composites  only  depends  upon  the  calibration 
of  the  bond  model  for  such  materials. 

The  research  was  carried  out  as  a  collaborative  effort  by  the  two  authors  and  constituted  the 
doctoral  research  of  the  first  author.  The  resulting  thesis  contains  additional  material  not  covered  in 
this  report,  such  as  a  description  of  the  finite  element  code  used  to  evaluate  the  numerical 
implementation  of  the  bond  model.  The  reader  is  referred  to  Mello  (1996)  for  this  additional 
information. 

The  research  had  four  main  components.  The  first  was  the  development  of  a  three- 
dimensional  composite  material  representation  of  a  reinforced  solid,  including  a  bond  model  that 
accounted  for  both  slippage  between  the  reinforcement  and  the  matrix  and  dilation  of  the  bond 
zone.  In  the  initial  development  the  bond  model  was  taken  to  be  linear  elastic.  The  composite 
material  representation  was  developed  by  using  a  micro-mechanics  model  consisting  of  an 
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approximate  "representative  volume"  for  an  assumed  hexagonal  packing  of  a  fiber-  or  bar- 
reinforced  solid. 

In  the  second  phase  the  resulting  composite  model  was  incorporated  into  a  plane  strain 
finite  element  code  with  four  degrees  of  freedom  at  each  node  (average  displacements,  bond 
slippage,  and  bond  dilation). 

The  third  phase  consisted  of  developing  and  calibrating  a  somewhat  simplified  version  of 
the  plasticity  bond  model  for  reinforced  concrete  that  was  originally  developed  by  Herrmaim  and 
Cox  (1992).  A  robust  numerical  implementation  scheme  for  evaluating  the  model  was  developed 
and  tested  for  a  range  of  bond  strain  and  stress  states.  The  bond  plasticity  model  was  then 
incorporated  into  the  composite  analysis  and  accompanying  two-dimensional  finite  element  code. 

The  final  phase  consisted  of  applying  the  analysis  to  a  number  of  tests  reported  in  the 
literature  for  reinforced  concrete.  The  results  of  this  final  phase  clearly  demcmstrated  the  feasibility 
of  incorporating  a  bond  model  that  includes  both  bond  slippage  and  dilation  into  a  composite 
representation  and  accompanying  finite  element  analysis  for  unidirectional  reinforced  materials. 

Structures  fabricated  of  reinforced  composite  materials  are  commonplace.  Reinforced- 
concrete  highway  bridges,  fiber-reinforced  polymer-matrix  aircraft  skins/spars,  and  high 
temperature  fiber-reinforced  ceramic-matrix  engine  components  serve  as  a  few  examples.  For 
modeling  purposes,  these  materials  can  be  viewed  as  consisting  of  three  general  constituents: 
fiber,  matrix,  and  bond.  Bond  refers  to  the  interaction  between  reinforcing  fibers,  or  bars,  and  the 
matrix  material.  For  many  reinforced  systems,  the  bond  behavior  is  very  important  to  the  overall 
response  of  the  structure. 

Because  these  composite  materials  are  used  in  many  important  structures,  the  need  to 
perform  detailed  structural  analyses  is  obvious.  Current  analyris  techniques  generally  fall  into  two 
categories:  composite  and  discrete.  In  a  composite  analysis,  the  reinforced  material  is 
approximated  as  a  continuum  with  properties  determined  fi’om  micro-mechanics  analyses  and/or 
empirical  test  data.  Mechanics  of  materials,  elasticity,  and  finite  element  analysis  techniques  for 
continuum  materials  may  then  be  applied.  These  analyses  generally  assume  perfect  bonding  and 
linear  response  of  fiber/matrix,  which  severely  limits  their  applicability.  They  may  be  valid  for  the 
analysis  of  fiber-reinforced  plastics  for  instance,  but  inadequate  for  reinforced  concrete  or  high 
temperature  ceramic-matrix  composites.  Discrete  analyses  employ  finite  element  techniques  to 
model  matrix,  fiber,  and  bond  separately.  Interface  elements  are  used  to  model  the  bond  between 
two-dimensional  or  three-dimensional  continuum  elements  which  separately  represent  the  matrix 
and  the  reinforcement.  If  nonlinear  material  models  exist  for  the  constituents,  then  the  results  of 
these  models  can  be  very  accurate.  The  problem  with  this  approach  is  cost;  it  is  costly  to  formulate 
and  develop  the  input  for  the  analyses,  and  computationally  expensive  to  run  the  very  large 
problems  that  result. 
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This  report  documents  the  development  of  a  composite  finite  element  that  includes 
nonlinear  response  due  to  the  interface  bond.  This  approach  results  in  finite  element  models  that 
are  relatively  easy  to  formulate,  inexpensive  to  run,  and  yet  capture  most  of  the  discrete  response 
behavior  of  the  constituents.  The  development,  which  includes  a  specialized  finite  element  scheme 
and  an  integrated  composite  micro-mechanics  analysis,  implements  an  elastic-plastic  bond  model 
that  was  developed  for  reinforced  concrete.  However,  the  resulting  analysis  technique  is  general 
and  may  be  applied  to  a  variety  of  reinforced  materials,  assuming  the  bond  model  used  in  the 
analysis  can  be  successfully  calibrated  for  the  system  at  hand. 


2.  REINFORCED  MATERIALS  AND  ANALYSIS 
2.1  Reinforced  Materials 

The  most  common  reinforced  material  is  a  unidirectional  composite.  This  composite 
configuration,  shown  in  Figure  1,  is  found,  for  example,  in  reinforced-concrete  bridge  structures 
and  aircraft  wing  spars.  A  reinforced  system  arranged  in  this  manner  forms  an  orthotropic 
composite. 

Microscopically,  this  material  is  heterogeneous  by  definition.  The  reinforcing  fibers  are 
typically  stronger  and  stiffer  than  the  matrix  material,  and  thus  the  resulting  longitudinal  stiffness 
and  strength  of  the  composite  is  usually  much  greater  than  in  the  transverse  directions.  The 
mechanical  properties  of  each  constituent  is  usually  well  characterized  and  documented. 

In  most  reinforced-material  structures,  the  load  is  applied  to  the  matrix  material  and  not 
directly  to  the  reinforcing  fibers.  The  load,  for  the  most  part,  is  transferred  to  the  fibers  through  a 
small  region  near  the  fiber  ends,  thus  creating  end  effects.  When  the  fiber  is  much  longer  than  the 
length  over  which  the  load  transfer  takes  place,  these  end  effects  are  usually  neglected.  In  the  case 
of  brittle  or  weak  matrix  materials,  such  as  ceramic-matrix  composites  or  reinforced  concrete, 
matrix  cracks  can  form,  and  these  end  effects  may  then  be  present  in  the  middle  of  the  structure. 

The  interfacial  bond  between  the  reinforcing  fiber  and  matrix  is  very  important  as  it 
influences  the  mechanical  properties  and  overall  performance  of  the  reinforced  material.  Because 
this  interface  is  responsible  for  transmitting  the  load  to  the  fibers,  the 
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importance  of  the  interface  has  been  recognized  and  much  research  has  been  undertaken  to 
understand  its  influence.  Most  researchers  consider  this  interface  to  be  really  a  third  material; 
Theocaris  (1987),  for  example,  calls  it  a  “mesophase”  and  has  devoted  an  entire  book  to  its  study. 

Because  of  the  complexity  of  incorporating  the  interface  behavior  into  the  analysis  of 
composite  materials,  the  bond  is  assumed  to  be  perfect  in  a  majority  of  cases,  and  designers  must 
often  rely  on  empirical  data  to  ensure  integrity. 

2.2  Reinforced  Materials  Analysis 

2.2.1  Perfect  Composite.  A  common  method  used  to  analyze  composite  materials  is 
to  assume  the  bond  between  the  fiber  and  matrix  is  perfect  and  that  no  local  movement  or 
deformations  exist.  If  one  further  assumes  that  the  reinforcement  is  arranged  uniformly,  and  that 
each  constituent  is  linear  elastic,  then  macroscopic  properties  can  be  obtained  by  analyzing  a  unit 
cell  or  representative  volume  element  (RVE)  of  material.  This  process  of  generating  composite 
continuum  properties  is  called  a  micro-mechanics  analysis.  Much  research  has  been  done  in  this 
area  (see  Hashin-Rosen,  1964  or  Jones,  1975,  for  example).  Most  micro-mechanics  analyses  use 
mechanics  of  materials,  elasticity,  or  energy  methods  to  compute  composite  elastic  constants  and 
strengths  in  terms  of  volume  fractions  of  constituents.  Aboudi  (1991)  in  his  text  on  micro¬ 
mechanics  uses  a  direct  approach  in  what  he  calls  the  "method  of  cells"  to  obtain  elastic  constants. 
He  has  also  examined  constituent  material  nonlineaiiiy  using  micro-mechanics.  This  work  will  be 
discussed  in  the  next  section.  Because  of  the  approximations  in  micro-mechanics  analyses, 
empirical  test  data  are  sometimes  used  to  supplement  the  theoretical  calculation  of  elastic  constants 
and  strength  data  (Tsai,  1992). 

Given  an  approximate  material  model,  derived  by  micro-mechanics  methods,  the  analyst 
may  use  well-established  mechanics  of  materials,  elasticity,  or  finite  element  techniques  to  analyze 
a  composite  structure  constructed  of  the  idealized  material.  This  method  works  well  for  polymer- 
matrix  composites,  such  as  graphite  fiber-reinforced  epoxy,  which  is  usually  limited  to  linear- 
elastic  deformations.  Care  still  must  be  taken  near  “close  outs”  or  “holes”  where  edge  effects  may 
be  present,  or  in  the  case  of  highly  loaded  structures  where  material  nonlinearities  are  indeed 
present 


2.2.2  Composite  Nonlinear  or  Imperfect  Interface  Micro-Mechanics. 
Aboudi  (1991)  has  developed  a  micro-mechanics  approach  to  model  composite  materials  with 
inelastic  matrix  materials  and  imperfect  bonding.  The  approach  is  an  extension  of  his  "method  of 
cells"  which  leads  to  average  composite  stress-strain  relations.  As  with  most  micro-mechanics 
methods,  the  reinforcement  is  assumed  to  be  in  a  periodic  array,  allowing  the  behavior  of  the 
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composite  to  be  found  from  analyzing  a  unit  cell  of  material.  Briefly,  his  analysis  of  the  unit  cell 
includes  imposition  of  traction  and  displacement  boundary  conditions  on  a  square  unit  cell.  The 
square  cell  is  divided  into  four  square  sub-cells  -  three  represent  matrix  material  and  one  represents 
the  fiber  material.  The  displacements  in  each  sub-cell  are  expanded  in  terms  of  distances  from  die 
center  of  the  unit  cell.  The  constituent  strains  are  evaluated  using  the  small  strain  tensor. 
Conditions  of  continuity  of  displacements  and  tractions  are  imposed  between  sub-cells,  and  an 
average  stress  theorem  is  used  to  derive  composite  stresses.  With  this  method,  Aboudi  has 
examined  sub-cell  interfacial  conditions  with  Mohr-Coulomb  frictional  behavior.  This  direct 
approach  has  advantages  in  that  it  can  also  be  extended  to  inelastic  composites,  while  a  minimum 
potential  energy-based  approach  cannot 

Aboudi  has  successfully  modeled  the  inelastic  behavior  of  metal-matrix  composites,  such 
as  graphite-aluminum  or  boron-aluminum,  for  loading  conditions  given  by  uniaxial  material- 
characterization  tests.  He  also  has  used  the  same  basic  approach  to  determine  effective  properties 
for  a  composite  whose  fiber-matrix  interfaces  are  completely  debonded.  These  general  methods 
are  considered  to  be  powerful,  yet  they  are  inadequate  in  modeling  an  interface  with  a  complicated 
plasticity  bond  model  that  depends  on  the  relative  axial  displacement  of  the  fiber  and  matrix  and  the 
dilation  of  the  bond  zone.  To  accommodate  such  complicated  interfacing  material  models  and  mote 
generalized  loading  conditions,  one  must  appeal  to  either  a  discrete  treatment  of  the  constituents  or 
a  hybrid  approach  that  treats  some  behavior  as  composite  while  handling  the  complex  material 
models  discretely.  A  discussion  of  both  follows. 

2.2.3  Discrete  Analyses.  For  reinforced-material  structures  where  edge  effects  or 
material  nonlinearities  are  important,  a  discrete  model  may  be  used  to  analyze  performance.  This 
method  employs  finite  element  techniques  to  model  separately  the  behavior  of  fiber,  matrix,  and  the 
bond  interface.  Continuum  elements  are  used  for  fiber  and  matrix,  and  an  interface  element  for  the 
bond.  See  Figure  2  for  a  sketch  of  what  a  finite  element  mesh  might  look  like  in  such  a  model. 

There  are  no  restrictions  with  this  type  of  analysis  except  those  associated  with  finite 
element  theory  itself.  The  matrix,  fiber,  and  bond  may  be  modeled  with  material  nonlinearities, 
assuming  that  adequate  material  models  exist.  The  problem  with  this  type  of  analysis  is  that  the 
cost  to  formulate  and  run  the  model  limits  the  applicability.  The  simple  reinforced  prism  in  Figure 
2  would  require  a  large  number  of  elements,  the  geometry  of  a  real  structure  would  be  significantly 
more  complex.  This  type  of  analysis  is  usually  reserved  for  a  sub-model  of  a  joint,  for  example, 
while  the  structural  designer  uses  a  simpler  more  approximate  analysis  for  the  complete  structure. 
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2.2.4  Composite  Discrete  Hybrid  Analysis.  There  have  been  some  analysis 
techniques  developed  that  lie  between  a  fully  homogenous  and  a  discrete  analysis  of  reinforced 
materials.  These  hybrid  analyses  offer  a  compromise  between  the  perfect  composite  assumption 
and  the  fully  discrete  treatment 

Pecknold  and  Hajali  (1993)  have  integrated  a  three-dimensional,  nonlinear,  micro¬ 
mechanics  materials  model  directly  into  a  structural  finite  element  analysis  for  laminated  composites 
using  standard  displacement-based  isoparametric  elements.  Their  analysis  technique  has  three 
levels:  (1)  a  nonlinear  micro-mechanics  material  model,  (2)  a  sub-laminate  model,  and  (3)  a  finite 
element  model.  The  resulting  twenty-node  brick  element  was  used  to  model  the  effects  of 
nonlinear  matrix  behavior  in  the  response  of  thick,  axisymmetrically  loaded  cylinders.  The  local 
stress  gradients  at  the  fiber-matrix  interfaces  were  not  treated  explicitly  with  this  procedure,  but  in 
other  respects  the  analysis  scheme  of  the  micro-mechanics  model,  integrated  with  the  finite  element 
procedure,  is  similar  to  that  used  in  the  present  study. 

Herrmann  and  Al-Yassin  (1978)  have  developed  a  composite  finite  element  analysis  for 
reinforced-soU  systems  whose  continuum  element  material  properties  reflect  the  response  of  the 
matrix  material  (soil),  the  reinforcing  sheets,  and  their  interaction.  This  composite  representation 
explicitly  models  reinforcement  slippage  between  constituents.  This  behavior  is  particularly 
important  for  soils  where  the  bond  results  firom  friction  and  is  weak.  The  technique  involved 
introducing  an  additional  nodal  degree  of  freedom  in  the  finite  element  scheme.  For  the  plane 
strain  case  that  was  studied,  the  unknowns  were  composite  x  and  y  displacements  and  slippage 
(the  relative  displacement  between  the  reinforcement  and  soil).  This  local  slippage  response  was 
approximated  as  a  continuous  field  variable.  This  approach,  combined  with  the  technique  of 
integrating  a  micro-mechanics  model  with  the  finite  element  scheme,  is  used  in  the  current  research 
involving  the  development  of  a  composite  element  that  includes  interface  bond  behavior. 


3.  COMPOSITE  FIBER.REINFORCED  ELEMENT  INCLUDING  BOND 
3.1  Objectives  and  Assumptions 

The  objective  of  this  work  was  to  develop  a  composite  finite  element  that  is  capable  of 
including  the  nonlinear,  elastic-plastic  response  of  the  bond  between  the  reinforcing  fibers  and 
matrix  material.  Such  an  element  is  applicable  to  reinforced-concrete,  ceramic,  or  carbon-matrix 
composites,  and  other  weak  or  brittle  matrix  materials  where  the  interface  with  the  reinforcement 
cannot  be  assumed  to  be  perfectly  bonded. 
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Several  important  assumptions  or  restrictions  ate  used  in  the  development: 

1.  Plane  strain  -  No  composite  strain  in  one  of  the  directions  transverse  to  the  fiber. 

2.  Unidirectional  reinforcement  -  Fibers  are  assumed  to  be  parallel  to  one  of  the  in-plane 
directions. 

3.  Linear  elastic  fiber-matrix  -  Fiber  and  matrix  constituents  are  assumed  to  be  isotropic 
and  linear  elastic. 

These  assumptions  allow  the  development  to  focus  on  inclusion  of  bond  without  the  added 
difficulties  of  three-dimensional  geometry,  angled  reinforcement,  or  nonlinear  fiber  and  matrix 
behavior.  These  effects  can  be  addressed  in  future  work  once  it  is  shown  that  the  elastic-plastic 
bond  response  can  be  included  in  a  composite  representation. 

3.2  Bond  Model  (Linear  Elastic) 

The  bond  model  that  is  the  key  feature  of  this  study  is  introduced  here.  It  is  a  simplified 
version  of  the  model  developed  by  Herrmann  and  Cox  (1992)  for  reinforced  concrete.  They 
studied  the  results  of  a  large  number  of  experiments  designed  to  determine  the  bond  behavior  of 
reinforced  concrete  and,  fi-om  this  information,  developed  a  comprehensive  bond  model.  The 
model  was  conceived,  tested,  and  shown  to  model  the  behavior  of  a  concrete-reinforcement 
interface.  The  linear-elastic  portion  of  the  bond  model  is  introduced  here,  as  it  is  essential  in 
understanding  the  composite  material  constitutive  relations  into  which  the  bond  model  is  to  be 
embedded.  The  details  of  the  simplified  nonlinear  model  are  described  later.  It  should  be  noted 
that  test  data  similar  to  that  used  to  develop  the  concrete  bond  model  exist  for  ceramic,  glass,  and 
metal-matrix  composites  (Makin,  Warren,  and  Evans,  1992;  and  Valente,  1994),  but  as  yet  no 
elastic-plastic  models  have  been  developed.  When  suitable  models  are  developed  for  these 
materials,  and  if  they  use  the  same  basic  framework  as  the  Herrmann  and  Cox  model,  the 
composite  finite  element  developed  here  could  then  be  used  to  model  structures  made  with  them. 

The  bond  model  idealizes  the  bond  zone,  the  interface  between  the  reinforcement  and 
matrix,  as  a  region  having  zero  thickness.  Aboudi  (1991)  made  similar  assumptions  when 
developing  micro-mechanics  relations  for  imperfect  interfaces.  The  bond  stresses  are  limited  to  a 
longitudinal  shear  stress,  x,  and  a  normal  stress  component,  a;  see  Figure  3.  The  conjugate  strains 
are: 
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qi  =  s/rb 

slip  strain 

(3.1) 

q2  =  6Al) 

dilation  strain 

(3.2) 

They  are  obtained  by  nonnalizing  the  axial  slip,  s,  that  is,  the  sliding  of  the  reinforcement  relative 
to  the  matrix,  and  the  dilation,  5,  of  the  bond  layer.  The  bar  radius,  rb,  is  used  to  render  the  bond 
model  independent  of  reinforcing  fiber  size. 

The  linear-elastic  portion  of  the  law  for  bond  developed  by  Herrmann  and  Cox  has  the 
following  form: 


t 
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(3.3) 


This  is  different  from  the  simpler  relation  used  by  Aboudi  which  ignored  dilation  of  the  bond  zone. 
The  complete  inelastic  bond  relations  derived  from  Herrmann  and  Cox’s  work  are  presented  in  a 
later  section. 


3.3  Composite  Element  Including  Bond 

The  composite  finite  element  scheme  includes  an  integrated  composite  micro-mechanics 
model.  These  composite  relations  are  developed  using  energy  methods  and  are  defined  in  detail  in 
Section  4.  The  resulting  three-dimensional  composite  linear  elastic  stress-strain  relations  are: 
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Bond  Shear  Stress 


Figure  3.  Bond  Model  Deformation  and  Stresses 


where  p  is  the  fiber  volume  fraction  and  the  stress-strain  matrix  [D]  is  a  function  of  Em,  Vm  (the 
matrix  elastic  constants),  Ef,  Vf  (the  fiber  elastic  constants),  [d]  (the  bond  elastic  constants),  and  p. 
The  definitions  of  the  stresses  and  conjugate  strains  are  developed  and  presented  in  Section  4 
following  the  definition  of  the  finite  element  displacement  unknowns. 

The  stress-strain  relations  given  above  are  used  in  a  finite  element  scheme  similar  to  that  of 
Herrmann  and  Al-Yassin  (1978)  discussed  earlier.  The  plane  strain  element  that  is  developed  is 
pictured  in  Figure  4.  The  displacement  unknowns  are  defined  as: 

u  -  X  displacement  of  matrix 

Uf  -  x  displacement  of  fiber 

V  -  composite  y  displacement 

6  -  bond  dilation 

Note  that  the  bond  slip  is  given  by 

S  =  U-Uf 

which  is  the  relative  displacement  of  the  fiber  and  matrix. 

Bilinear  shape  functions  are  used  to  describe  both  matrix  and  fiber  displacement  as  nodal 
unknowns.  Fiber  displacement  (instead  of  slip)  is  selected  as  a  global  unknown  because  it  offers 
more  flexibility  in  imposition  of  boundary  conditions. 

The  slip  and  dilation,  which  are  local  to  the  interface  of  fiber  and  matrix,  are  modeled  as 
continuous  field  variables  as  in  the  work  of  Herrmann  and  Al-Yassin  (1978).  With  fiber 
displacement  or  slip  as  an  unknown,  one  can  account  for  bond  breakdown  and  still  ensure 
compatible  inter-element  fiber  displacements.  The  use  of  dilation  as  a  field  variable  is  less 
straightforward;  it  is  selected  as  an  unknown  as  it  simplifies  implementation  of  the  nonlinear  bond 
model.  It  is  recognized  that  for  the  linear  elastic  case  it  could  be  condensed  out  at  the  element  level 
if  one  were  willing  to  sacrifice  dilation  inter-element  compatibility.  This  may  not  be  a  bad 
assumption  in  light  of  other  approximations  made  in  the  micro-mechanics  material  model. 

Given  the  composite  constitutive  relations  and  nodal  unknown  definitions,  one  can  use 
well-established  finite  element  formulation  techniques  such  as  the  principle  of  virtual  work  or 
minimization  of  potential  enCTgy  to  develop  a  composite  continuum  finite  element  It  is  also  shown 
that  this  work  can  be  extended  to  the  nonlinear  case  of  an  elastic-plastic  bond  model  using  the  same 
basic  framewoik.  Details  of  this  development  are  discussed. 
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Figure  4.  -  Composite  Plane  Strain  Global  Unknowns 


4.  DEVELOPMENT  OF  LINEAR-EQUIVALENT  HOMOGENEOUS 
COMPOSITE  PROPERTIES  INCLUDING  BOND 

The  finite  element  scheme  requires  an  equivalent  homogeneous  material  model  as  discussed 
in  Section  3.  The  basic  goal  is  to  find  effective  elastic  moduli  of  the  reinforced  material  in  terms  of 
the  constituents  elastic  moduli  and  their  geometry.  The  details  of  the  development  of  a  composite 
material  model  that  consists  of  linear-elastic  fiber,  matrix,  and  bond  follow. 

Aboudi  (1991),  Bienveniste  (1985),  and  others  have  performed  micro-mechanics  analyses 
with  imperfect  interfaces  or  elastic  tangential  slip.  The  micro-mechanics  analysis  of  this  study  is 
different  in  that  it  includes  the  bond  dilation  effect  in  addition  to  slip. 

The  general  method  used  is  an  extension  of  Hashin  and  Rosen's  (1964)  work. 
Deformation  states  of  a  representative  volume  element  of  material,  whose  strain  energies  are 
tmcoupled,  are  examined.  In  each  case,  homogeneous  properties  are  obtained  by  equating  the 
virtual  work  of  the  composite  system  to  that  of  the  constituents.  The  virtual  work  approach  is 
selected  over  equating  strain  energies  because  it  is  more  general,  allowing  for  a  nonsymmetrical 
bond  model,  and  it  is  also  easily  extended  to  the  case  of  a  nonlinear  inelastic  bond  response,  which 
is  discussed  later. 

4.1  Hexagonal  Pack  Representative  Volume  Element 

The  fiber  and  matrix  are  assumed  to  be  arranged  in  hexagonal  pack  as  shown  in  Figure  5. 
This  packing  assumption  was  introduced  by  Hashin  and  Rosen.  It  allows  one  to  consider  the 
response  of  a  single  cylindrical  representative  volume  element  (RVE),  which  is  a  sample  of  the 
composite  material  that  represents  the  composite  on  average.  It  contains  the  three  phases  of 
interest:  fiber,  bond,  and  matrix.  The  RVE  is  large  compared  to  the  constituent  micro-structure, 
but  small  relative  to  the  entire  structure.  This  allows  one  to  use  continuum  theories  for  scales 
larger  than  the  RVE.  Because  the  fiber  lengths  are  much  greater  than  the  fiber  diameters,  the 
analysis  can  be  further  simplified  by  examining  a  very  long  cylindrical  sample  of  matp.riaf 

At  this  point,  the  response  of  the  fiber,  matrix,  and  bond  is  assumed  to  be  linear  elastic, 
thus  the  constituent  material  properties  for  the  micro-mechanics  model  are: 

Em»Vin  -  Young's  modulus  and  Poisson's  ratio  for  matrix  material 

Ef,  Vf  -  Young's  modulus  and  Poisson's  ratio  for  reinforcing  fiber  material 
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[D]  -  stress-strain  matrix  for  the  composite  with  elastic  bond  behavior 


The  composite  section  geometry  parameter  is  defined  using  the  fiber  volume  fraction 

_ 

P  7 

Jta^ 


(4.1) 


for  the  representative  volume  element,  where  is  the  fiber  radius  and  "a"  is  the  outer  diameter  of 
the  RVE.  The  ratio  of  the  RVE  outer  radius  to  fiber  radius  is  also  used  in  the  development  and  is 
defined  as 


4.2  Unit  Cell  Deformations,  Strain  and  Stress  Definitions 


(4.2) 


Because  of  the  peculiarities  of  the  analysis  that  include  bond  slip  and  dilation,  it  is 
important  to  carefully  define  the  RVE  deformations  that  are  shown  in  Figure  6.  The  displacements 
are: 


u,v  -  average  transverse  displacement 
w  -  average  matrix  displacement  in  the  fiber  direction 
uf  -  fiber  displacement 
s  -  slip 

5  -  average  dilation  of  the  bond  zone 
The  composite  strains  that  result  from  the  micro-mechanics  analysis  are: 
Sxj  ey»  Yxy*  Yxz,  Yyz  '  composite  strains 

Ez  -  average  matrix  strain  (fiber  direction) 

Ef  -  average  axial  fiber  strain 
qi  -  average  slip  strain 
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Figure  6.  -  RVE,  Unit  Cell  Deformations 
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The  conjugate  stresses  are: 
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-  composite  stresses 

average  matrix  stress  (fiber  direction) 

-  average  axial  fiber  stress 

-  average  bond  zone  shear  stress 

-  average  bond  zone  normal  stress 

-  fiber-bond-matrix  equilibrium  stress 


The  average  stresses  and  strains  are  for  one  constituent,  but  are  averaged  over  the  entire  unit  cell. 
The  composite  properties  are  also  average  values,  but  are  representative  of  the  combination  of 
constituents.  Six  deformation  cases  are  used  to  determine  the  effective  elastic  constants  in  the 
matrix  [D].  In  each  of  these,  a  deformation  is  applied  that  results  in  a  strain  energy  that  is 
uncoupled.  These  deformations  are  used  to  determine  boundary  conditions  for  the  field  problems 
to  be  solved  for  the  RVE.  The  solutions  of  these  problems  are  used  when  equating  the  virtual 
work  of  the  homogeneous  material  to  that  of  the  constituents. 

4.3  Case  A  -  XZ  Shear 


The  first  deformation  case  examined  yields  qo,  ait  effective  shear  modulus  for  pure  shear  of 
the  RVE  in  the  xz  plane.  Since  no  strain  energy  is  developed  at  the  bond  surface,  bond  stress  and 
strain  does  not  enter  this  case.  This  deformation  is  shown  in  Figure  7.  For  the  composite 


'Cxz  =  qoYx2  (4.3) 

where  Xxz  and  are  composite  stress  and  strain;  the  variable  qo  is  the  effective  shear  modulus 
sought. 
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The  prescribed  deformation  is  Yxz’  while  all  other  strains  are  zero  for  this  case.  The  axial 
displacement,  Uz ,  at  the  boundary  is  approximated  as: 

Uzlr=a  =  f(a)  cos(0  )  (4.4) 

and  similarly  the  fiber  and  matrix  z  direction  displacements  are  written  as: 

Uzf  =  ff(r)  cos(0 )  (4.5) 

Uzm  =fm(r)cos(0)  (4.6) 

The  subscript  "f '  denotes  fiber  and  "m"  matrix.  The  above  boundary  condition.  Equation  4.4,  is 
applied  to  Uzm  which  results  in : 

fin(a)  —  u  Yxz  (4-7) 

FibCT  and  matrix  interface  compatibility  and  equilibriinn  require: 

ff(rb)  =  fni(rb)  (4-8) 

'^xzjji  (tb)  ~  '^xzf  (ti))  (4.9) 

The  deformation  for  either  fiber  or  matrix,  Equations  4.5  and  4.6,  is  written  as: 

Uzi  =fi(r)cos(0)  (4.10) 

where  i  is  either  m  or  f.  The  strain-displacement  relations  in  cylindrical  coordinates  reduce  to: 

£rrj  =  £69i  =£zzj  =Yrzj  =  0  (4.11) 

Yrzi=^cos0  (4.12) 

Y0zi=^sin0  (4.13) 
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Hooke’s  law  for  the  fiber  and  matrix  materials  is  written  as: 


II 

o 

(4.14) 

'^ezj  =  Gi  Yezj 

(4.15) 

The  equilibrium  equation  reduces  to: 

^  ^  1  ^ezi  1 
dr  r  90  r 

(4.16) 

Equations  4.12,  4.13,  4.14,  4.15,  and  4.16  form  a  complete  set  of  elasticity  equations.  The 
solution  to  these  equations  is  of  the  form 

fi  =  rc  (4.17) 

and  is  subject  to  the  boundary  conditions  of  Equations  4.7, 4.8,  and  4.9,  and  the  further  condition 
that  the  result  is  finite  at  r  =  0.  The  displacement  solutions  are: 


Uzf  =  Af  rcos0 

(4.18) 

Uzm  =  (Amr  +  ^)cos0 

(4.19) 

where  the  arbitrary  constants  Af,  Am,  and  Bm  are  found  firom  the  given  boundary  conditions.  The 
results  are  found  to  be: 


Am  =  KiYx2 

(4.20) 

Bm  =  a2K2Yx2 

(4.21) 

Af  =  K3Yxz 

(4.22) 
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where 


(Gm/Gf+l)a2 

^  (a2-l)+(l+a2)Gm/Gf 

(4.23) 

K> 

II 

1 

(4.24) 

K3  =  [a2(l-Ki)  +  Ki] 

(4.25) 

with  a  defined  previously  in  Equation  4.2  as  the  ratio  of  the  matrix  radius  to  fiber  radius.  Gm  and 
Gf  are  the  elastic  shear  moduli  for  the  matrix  and  fiber,  respectively.  Given  the  displacement 
solution,  one  can  determine  strains  and  stresses  for  each  constituent  using  Equations  4.12, 4.13, 
4.14,  and  4.15.  To  derive  an  expression  for  the  desired  effective  homogeneous  shear  modulus, 
the  virtual  work  of  the  equivalent  composite  homogeneous  material  is  equated  to  that  of  the 
constituent  materials,  i.e.. 

Wc  =  Wm  +  Wf 

(4.26) 

where  Wc,  the  composite  virtual  work,  is  defined  as: 

Wc  =  Xxz  Ayx2  rdr  d0 

(4.27) 

and  where  Ay^z  is  a  virtual  strain  applied  to  the  composite  RVE.  Likewise  for  the  matrix  and 

fiber: 

Wm=  [tr^AY,z„  +  Te2„ATe2„]rdr<ie 

(4.28) 

Wf  =  [Tjzf  Ayrzf  +  XBzf  ATBzflr  dr  dO 

(4.29) 

The  AVrzf,  Ayezf,  AYrzjjj,  and  Ayez^jj  virtual  strains  are  written  in  terms  of  Ayx^.  It  has  also  been 
shown  that  Xizp  XGzp  ^nd  xqz^  can  be  written  in  terms  of  yxz  and  the  elastic  constants  from 


the  elasticity  solution.  Equating  virtual  work,  as  described  in  Equation  4.26,  one  can  obtain  the 
desired  composite  shear  modulus,  "no: 

TlO  =  +  Gm  [Ki2  (1  -  -i- )  +  K22(a2  - 1)]  (4.30) 

a2  a2 

Note  that  this  homogeneous  or  composite  shear  modulus  is  a  function  of  the  elastic  constants  of  the 
constituents  and  the  volume  fraction  of  each  in  the  RVE. 

4.4  Case  B  -  YZ  Shear 

The  next  deformation  case  examined  is  pure  shear  in  the  yz  plane.  Because  of  the 
axisymmetry  of  the  RVE,  the  results  of  this  case  are  the  same  as  for  Case  A 


'Cyz  =  'noYyz  (4.31) 

where  Tio  is  given  by  equation  Equation  4.30.  Bond  stress  and  slip  are  absent  from  this  case  also. 

4.5  Case  C  -  Change  in  Volume  with  Slip  and  Dilation 

The  next  case  considers  the  following  deformations: 

s  -  slip 

S  -  dilation 

Ef  -  fiber  axial  strain 

Ez  -  matrix  strain  (in  fiber  direction) 

Ex  -  transverse  strain 
Ey  -  transverse  strain 

where  Ex  =  £y  =  ei/2,  which  corresponds  to  a  change  in  diameter  of  the  RVE.  The  deformation 
state  is  pictured  in  Figure  8.  Recall  how  the  strain  for  homogeneous  isotropic  materials  may  be 
decomposed  into  hydrostatic  and  deviatoric  components.  The  present  deformation  case  is 
analogous  to  the  hydrostatic,  or  change  in  volume,  component 


23 


The  inclusion  of  dilation  as  a  global  degree  of  freedom  requires  that  one  of  the  constituent 
relations  developed  in  this  deformation  case  be  a  statement  of  equilibrium  at  the  interface.  A 
fictitious  interface  stress,  a*,  is  introduced  conjugate  to  the  dilation  strain.  To  aid  in  understanding 
this  procedure,  a  simple  analogous  one-dimensional  example  is  reviewed  prior  to  introducing  the 
much  more  complicated  case  (see  Figure  9).  For  this  simple  bar  problem,  the  dilation  of  the 
interface  at  the  top  of  the  bar,  5,  and  the  deflection  of  the  end  of  the  bar,  u,  are  the  displacement 
unknowns.  Virtual  work  is  used  to  find  the  constitutive  relations  for  F*  and  F,  the  conjugate 
stress  quantities.  If  one  then  studies  F*,  it  is  seen  that  it  is  an  equilibrium  expression  for  the 
interface. 

The  bar  strain  and  virtual  strain  are  written  as: 


e  =  (u-8)// 

(4.32) 

Ae  =  (Au  -  AS)// 

(4.33) 

For  simplicity,  the  bar’s  area  is  assumed  to  be  unity  and  the  interface  constitutive  relation  is  taken 
to  be: 


Fb  =  Kb  S  (4.34) 

The  composite  virtual  work  is  equated  to  that  of  the  constituents  as  is  in  the  more  complicated  case, 

Wc  =Wbar  +  Wbond 
Wc  =FAu  +  F*  AS 
Wbar  =J^  AeEedx 

Wbond  =  Kb  8  (AS) 

where  E  is  Young’s  modulus  for  the  bar. 
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Figure  9.  -  Simple  Bar  Interface  Example 


Combining  the  above  expressions,  simplifying  the  results,  and  substituting  the  appropriate 
strain  expressions,  one  finds: 

[F  -  y  (u-5)]Au  +  [F*  +  j  (u-5)  -  Kb  6]A5  =  0 

For  arbitrary  Au  and  AS,  it  follows  that: 

F=  I  (u-5) 

which  is  the  constitutive  relation  for  the  bar  portion  of  the  composite  in  terms  of  u  and  5,  and 

F*=E^-Kb6  (4.37) 

which  upon  inspection  reveals  that  F*  should  be  zero  for  equilibrium  at  the  interface.  One  may 
question  where  to  set  F*  to  zero  to  ensure  equilibrium  if  the  constituent  relations  of  Equations  4.36 
and  4.37  were  to  be  used  in  a  finite  element  scheme.  To  show  how  this  is  achieved,  one  simply 
writes  the  expression  for  potential  energy  for  the  bar-interface  composite,  and  the  resulting  energy 
mhumization  process  accompanying  the  finite  element  development  results  in  F*  being  set  to  zero. 
(Recall  that  the  finite  element  equations  are  derived  by  taking  derivatives  of  the  potential  energy 
with  respect  to  u  and  5.) 

Having  introduced  the  concept  of  a  fictitious  stress  at  the  interface,  we  may  now  proceed 
with  analyzing  the  deformation  case  pictured  in  Figure  8.  The  prescribed  RVE  deformations  arc  s, 
5, ef,ez  and 


(4.35) 


(4.36) 


Ex  =  Ey  =  ei/2  (4.38) 

The  fiber  displacement  is  uf  =  w  -  s,  and  the  fiber  strain  is: 

ef=e2-^  (4.39) 

In  the  radial  displacement  of  the  boundary  of  the  RVE,  d  needs  to  be  related  to  toe  prescribed  strain 
in  the  x  and  y  directions  given  above.  To  calculate  d,  one  equates  the  changes  in  volume: 
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%  s?  (Ex  +  Ey)  =  27ca(d) 


(4.40) 


Using  Equation  4.38,  one  finds 


d 


(4.41) 


which  is  the  desired  boundary  displacement  in  terms  of  the  radial  strain. 

As  in  Case  A,  one  uses  the  elasticity  relations  for  a  long  axisymmetric  cylinder  to  calculate 
each  constituent’s  stress  and  strain  that  result  firom  the  prescribed  RVE  deformations. 

Conditions  of  axisymmetric  deformation  and  the  RVE  being  long  lead  to  the  assumptions 
of  no  changes  in  the  6  and  z  direction.  Radial  direction  equilibrium  then  gives: 


dOr  .  Or-ae 
- 


=  0 


(4.42) 


The  strain  displacement  relations  (with  constant  e^  )  are: 


Er  = 


du 

dr 


ee  = 


Ik 

r 


(4.43) 


Hooke’s  Law  for  the  assumed  linear-elastic,  isotropic  material  is: 
(Tr  =  (X  +  2|l)Er  +  X  E0  +  XEz 
00  =  (X  +  2p)E0  +  X  Er  +  XEz 


(4.44) 


The  displacement  formulation  is  obtained  by  substituting  Equations  4.43  into  4.44  and  then  the 
subsequent  results  into  the  equilibrium  Equation  4.42 


d^Ur  .  i  ^  it  _n 
dr2  r  dr  ■  1-2  "  ” 


(4.45) 


which  holds  for  either  matrix,  m,  or  fiber,  f.  This  second  order,  ordinary  differential  equation  has 
a  general  solution  of  the  form: 
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Ur  =  cir  +  ^ 


(4.46) 


The  radial  displacement  for  the  matrix  is: 

The  radial  displacement  for  the  fiber  reduces  to: 

Urf  =  cif  r  (4.48) 

because  the  solution  must  be  bounded  at  r  =  0.  The  strain  displacement  relations  are  written  in 
terms  of  these  displacement  functions  and  the  prescribed  axial  deformations.  For  the  matrix: 

=  (4.49) 

ezm  =  ez 
and  for  the  fiber 

—  Clj 

e0f  =  cif  (4.50) 

ds 

ezf-ef=ez-  ^ 

The  boundary  and  interface  conditions  are  now  defined  to  complete  the  elasticity  problem.  For  the 
boundary  of  the  RVE  firom  Equation  4.41: 

«rmla  =  ‘i  =  af  (4-51) 
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For  the  interface,  equilibrium  and  continuity  yield: 


where  5  is  the  average  RVE  bond  dilation. 

This  elasticity  problem  is  now  solved  for  cij„,  C2ni>  and  cif.  Mathematica  (1995)  was 
employed  to  check  the  algebra  involved  in  the  solutions  with  the  following  results: 


cif  =  Keei  +  K75/rb  +  KgCz  +  K9ef 
cim  =  Kioei  +  Kii5/rb  +  Ki2ez  +  Kisef 
C2tn  =  a2(Ki4ei  -  Kii5/rb  -  Ki2ez  -  Ki3ef) 


where 


K4  =  (Ma2+l)  +  Xm)/(a2-l) 

Ks  =  l/(2(K4  +  lif+?if)) 

K6  =  a2(K4-llm)K5 

K7  =  -2K4K5 

Kg  =XmK5 

K9  =  -^fKs 


Kio 


a2-2K6 

2(a2-l) 


(4.53) 


(4.54) 


Kii 


-(I+K7) 

(a2-l) 


Ki2  =  -K8/(a2-l) 
Ki3  =  -K9/(a2-l) 
Ki4  =  2  -KiO 
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The  algebra  makes  these  expressions  appear  complicated  but  once  again  they  are  simply  functions 
of  the  constituent  elasticity  constants  X  and  p,  and  the  quantity  a  which  is  related  to  the  volume 
fraction  of  each  of  the  constituents. 

What  remains  is  to  equate  the  virtual  work  of  the  desired  composite  model  to  that  of  the 
constituents.  For  this  case: 


Wcomp  —  Wfiber  +  Wmatrix  +  Wbond 


(4.55) 


For  the  composite,  the  work  is: 

Wcomp  ~  ^3.^  { <5yAey  + 

(l-p)OzAe2  +  pOfAef  +  (x  As  +  c*  A6)}  (4.56) 

Tca^ 

Some  important  things  to  note  in  the  preceding  equation  are  the  volume  fraction  quantities  in  front 
of  the  matrix,  fiber,  and  bond  terms  and  that  the  x  and  y  subscripted  terms  are  for  the  composite. 
Note  the  inclusion  of  a*,  the  fictitious  composite  bond  stress,  introduced  previously  with  the  help 
of  the  simple  bar  example.  Recall  that  6x  =  £y  =  ei/2,  and  so  Equation  4.56  now  is  written  as: 

Wcomp  =  7ia2{oi Aei  +  (l-p)az  AEz  +  pGf  A£f  +  (x  As  +  a*  A5) }  (4.57) 

7ca2 


The  fiber  virtual  work  is  written  as: 


Wfiber  =  (<yrfA£rf  +  O0fA£ef  +  OfAef)  rdidG 


(4.58) 


The  bond  virtual  work  is: 


Wbond  =  liUTb  (x  As  +  aA5) 


(4.59) 


The  matrix  virtual  work  is: 
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Wmatrix  =  «ird0 


(4.60) 


TTie  constituent  virtual  strains  in  these  expressions  are  found  using  the  strain  relations  of  Equations 
4.49  and  4.50.  For  example,  for  the  radial  strain  in  the  matrix,  one  finds: 

=  Aciixi  -  Ac2m  ^ 

A5 

Acii^  =  Kio  Aei  +  Kii  —  +  K12  Aez  +  K13  Acf 

Ac2„  =  a2(Ki4Aei-Kii  ^  -  K12 AEz  -  K13  Aef) 

and  so  on  and  so  forth  for  Aeg^^,  ACr^,  and  Ae0j. 

Mathematica  was  once  again  employed  to  do  the  algebra  and  find  the  resulting  constitutive 
relations;  details  are  to  be  found  in  Mello  (1996).  The  final  results  are: 


ai=Tiiei+Ti2ez+Tl3ef+Tl4  8/rb  (4.61) 

Oz  =  [1)2  £1  +  T15  fz  + 116  Ef  +  717  S/rbl  (4.62) 

(1-p) 

Of  =  -  [113  El  +  716  Ez  +  718  Ef  +  TI9  S/n,]  (4.64) 

p 

t  =  [diis/rb  +  di2  6/rb]  (4.65) 

<7*  =  •%  [714  El  + 117  Ez  +  718  Ef  +  TllQ  5/rb] 
a2 

+  [d2i  s/rb  +  d22  5/rb]  (4.66) 


As  explained  before,  the  quantity  c*  is  an  equilibrium  equation  for  the  interface  and  is  set  to  zero 
during  the  process  of  developing  the  finite  element  equations  for  the  composite  element  based  on 
these  material  relations.  The  Tj's  are  the  desired  effective  moduli;  detailed  expressions  are  to  be 
found  in  the  Mathematica  Notebook  given  in  Mello  (1996).  (The  expressions  are  in  FORTRAN 
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format  and  were  pasted  directly  from  Mathematica  into  the  FORTRAN  finite  element  program.) 
The  e)q)ression  for  T)  i  given  below  is  representative  of  these  results: 

T]!  =  4  (jinri  K^q  -  [Xm  +  ^.m 

^  M'in  frjj  -  r^j  Pin  KjQ/a^ 

+  r^pfK^/a2  +  rJXfK^a2 

-r2XmKj/a2 


The  parameter  Tii  given  above,  as  are  all  the  others,  is  a  function  of  Pm»  ^m.  W,  ib  and  "a" 
only. 

Finally,  note  that  Equation  4.61  relates  the  composite  stress,  oi,  to  the  composite 
deformations.  Recall  that  ai  and  ei  are  the  radial  stress  and  strain.  What  is  desired  ultimately  are 
Ox  and  Oy.  These  are  found  by  superimposing  the  results  of  Cases  C  and  D,  which  follow. 

4.6  Case  D  >  Cross  Section  Distortion 

The  next  deformation  case  is  shown  in  Figure  10;  for  this  case  the  fiber  strain,  ef ,  and  the 
matrix  strain  (in  the  fiber  direction),  Ez  ,  are  zero,  which  also  implies  no  sUp.  The  dilation  is  also 
taken  to  be  zero.  Recall  that  this  is  the  average  value  for  the  RVE.  The  cross  section  is  distorted 
such  that 


ex  =  -ey  =  e2/2  (4.67) 

where  £2  is  the  imposed  strain,  and  the  deflections  in  the  x  and  y  directions  are: 

£2 

Ux=y  X 

Uy  =  -Yy  (4.68) 
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Y 


Figure  10.  -  Case  D,  Cross  Section  Distortion 


This  deformation  is  similar  to  the  direct  components  of  deviatoric  strain  when  a  material  response 
is  divided  into  hydrostatic  and  deviatoric  components.  Recall  that  Case  C  was  analogous  to  a 
change  in  volume  strain.  This  deformation  of  the  RVE  is  the  nonaxisymmetric  plane  strain 
behavior  of  a  thick  cylinder.  In  terms  of  cylindrical  coordinates,  it  is  observed  that  the  radial 
displacement,  %,  is  a  symmetrical  function  of  0  and  the  tangential  displacement,  ue,  is  an 
antisymmetric  function  of  0.  For  this  type  of  periodic  behavior  one  seeks  solutions  of  the  form 

%  =  f(r)  cos(n0) 

ue  =  g(r)  sin(n0)  (0  <  0  <  2k)  (4.69) 

where  n  =  2.  The  deflection  at  the  boundary  of  the  RVE  likewise  is  approximated  as: 


Ur|^  =f(a)cos(20) 
ue|  =g(a)sin(20) 


(4.70) 


The  Ur  and  ue  displacements  are  expressed  in  terms  of  the  x  and  y  coordinates  and  the  prescribed 
deformations  of  Equation  4.68  with  the  appropriate  transformation  to  cylindrical  coordinates, 
which  is: 


x=  rcos(0) 
y=  rsin(0) 


(4.71) 


The  radial  displacement  is  written  as: 


e2 


Ur=  rcos(20) 


(4.72) 


using  the  appropriate  substitutions  and  trigonometric  identities.  For  the  outer  boundary  of  the  RVE 
at  r  =  a,  one  finds 


Ur  I  =  y  acos(20) 


(4.73) 


which  implies  that 


f(a)=f  a 


(4.74) 


on  the  boundary.  One  can  likewise  show  that: 


U0  =  -  Y  rsm20 


(4.75) 


and 


g(a)  =  -  y  a 


(4.76) 


The  equilibrium  equations  for  fiber  and  matrix  for  the  elasticity  problem  associated  with 
this  deformation  case  are: 
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(4.77) 


^  r  30  r 


^Trt=o.o 

dr  r  30  r  ™ 


(4.78) 


Likewise  the  strain-displacement  relations  are: 


dur 

^r=-3F 


ee=i^  +  St 

^  T  r 


r  30  +  ar  r 


ez  =  Y0z  =  Yrz  =  O 


(4.79) 


for  both  the  fiber  and  the  matrix.  The  constitutive  relations  are: 


. . [(i-v)er+vee] 

(l+v)(l-2v) 


^6=  ,,  ,,,  ^  ,  [(l-v)ee-hver] 

(l+v)(l-2v) 


=  ,,  ,,,  ,  [V  (er  +  ee)] 

(l+v)(l-2v) 


'Z'. — : 
2(l+v) 


X0Z  —  Xiz  —  0 


(4.80) 


where  E  and  v  are  the  appropriate  Young's  modulus  and  Poisson’s  ratio  for  fiber  or  matrix 
depending  on  the  solution  sought  The  displacement  solution  for  the  matrix  is  found  to  be: 


fm  =  Cl  r3  +  C2„  r  +  C3„  7  +  C4„.  4 

*“  4|n  .i-m  ->ni  j-  tiji  ^.3 
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(4.81) 


gn.  =  -(^)Cl.r3.C2,r.  C3„  i  .C44 


2(l-v) 


and  for  the  fiber  whose  solution  must  be  bounded  at  r=0: 


ff  =  Cifr3+C2fr 


,  3-2v  ^  ^  ^ 

gf=-(— —  )Cifr3-C2fr 

2v  ^  ^ 


(4.82) 


The  boundary  condition  at  the  RVE  outer  radius  "a"  is  from  Equations  4.74  and  4.76: 
fm(^)  “  ^ 

gm(a)  =  -ay  (4.83) 

Assuming  perfect  bond,  the  interface  conditions  are: 

Urf(rb)  =  Ur^(rb) 

U0f(rb)  =  u0jjj(rb) 

(rb)  =  (^b) 

'Cr0f(rb)  =  Tr0^(rb)  (4.84) 

These  conditions  result  in  six  equations  which  are  used  to  solve  for  the  constants  Ci^^, 

Qm’  ^2f  of  the  displacement  solution  given  in  Equations  4.81  and  4.82. 

Attempts  to  solve  this  set  algebraically  with  tire  symbolic  program  Mathematica  resulted  in 
terribly  complicated  algebra.  Since  all  results  are  to  be  implemented  numerically,  in  a  FORTRAN 
program,  it  was  decided  to  leave  the  solution  of  this  system  to  a  munerical  linear  algebra  technique 
such  as  Gaussian  elimination. 

As  in  the  other  cases  previously  studied,  the  effective  elastic  moduli  are  obtained  by 
equating  the  composite  virtual  work  to  that  of  the  constiments  for  the  prescribed  deformation: 
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Wc  =  Wf  +  Wm 


(4.85) 


For  the  composite 

”  7t3.^  [^x  ^^x 

or 

Wc  =  7ia2a2^  (4.86) 

for  the  deformation  prescribed,  where  A  again  denotes  virtual  quantities.  For  the  fiber  the  virtual 
work  is: 

=  [Off  Aetf +  oef  Aeof+Tref  A7r0f]  rdrdG  (4.87) 

and  for  the  matrix: 

Wm  =  [<Jr„  +  <50^  r,0^)  rdrde  (4.88) 

Equating  the  virtual  work  is  carried  out  in  terms  of  the  six  arbitrary  constants  (which  are 
determined  numerically).  The  result  is: 

<y2  =  'nilE2  (4.89) 

where  1)11  is  the  effective  modulus  sought  for  this  deformation  case.  The  modulus  ilii  is  found  to 
be: 

^11=  (PlCl„2  +  P2Ci„C2,+  |i3Ci„C2, 

+  P4Cif2  +  p5CifC2f+P6C2f2 

+  P7C3,C4,  +  P8C4,2.,P9C3„2 

+  PioC2^2)/d  (4.90) 
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where  D  is 


D  =  2a8  Vf2  (1  +  Vf)(Vm  -  1)2  Vn,2  (1  +  Vr,)  (4.91) 

and  the  betas  are 

Pi  =  3a6  (-1  +  a6)Ein  Vf2  (1  +  Vf)(3  -  2  Vm)(-1  +  Vni)2 
p2  =  12a6  (-1  +  a4)Em  Vf2  (1  +  Vf)(-1  +  Vm)2 
Ps  =  6a6  (-1  +  a2)Ein  Vf2  (1  +  Vf)  Vm  (1  -  3Vm  +  2Vm2) 
p4  =  3a6  Ef  (3  -  2Vf)(-l  +  Vm)2  Vm2  (1  +  Vm) 

Ps  =  12a6  Ef  Vf  (- 1  +  Vm)2  Vm2  (1  +  Vm) 

P6  =  8a6  Ef  Vf2  (-1  +  Vm)2  Vm2  (1  +  Vm) 

P7  =  12a2  (-1  +  a4)Em  Vf2  (1  +  Vf)(l  -  Vm)  Vm2 
p8  =  24  (-1  +  a6)Em  V^  (1  +  Vf)(-1  +  Vm)2  Vm2 
p9  =  a4  (-1  +  a2)Em  Vf2  (1  +  Vf)(3  -  2  Vm)  Vm2 

PlO  =  8a6(-l  +  a2)Em  Vf2(l  +  Vf)(-1  +  Vm)2  Vm2  (4.92) 

which  are  basically  the  coefficients  of  the  arbitrary  constant  C’s.  The  six  equations  that  result  from 
Equations  4.81  and  4.82  are: 

Cif  +  C2f  =  +  C2^  +  CSjjj  +  (4.93a) 
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Ci£  -  C2f  -  3  Cij^  /  2vf  = 


-  C2jjj  +  +  Cijjj  (-3  +  2  Vin)/2Vni 

+  C3„  (-1  +  2Vm)  /  (2  -  2  Vm)  (4.93b) 

C2fEf/(l  +Vf)  = 

Em  (-C2^  +  C3„  +  3C4^  +  C2^  Vm  -  3  C4^  Vm)  /  (-1  +  Vm^)  (4.93c) 
-Ef  (3  Cif  +  2C2f  Vf)  /  (2  Vf  (1  +  Vf))  = 

Em  ('3  Cijjj  +  3  Cijjj  Vm  -  2  C2j^  Vm  -  03^^^  Vm 

-  6  C4„  Vm  +  2  C2^  Vm2  +  6  €4^  Vm^)  /  (2  Vm  -  2  Vm^)  (4.93d) 

®  ^2m  +  +  Tb^  C3jn/a  +  rb4  C4ja?  =  a  E2/2  (4.93e) 

■a  C2jjj  +  rb4  Cl4jjj/a^  +  a^  ^  Vm) 

+  C3„  11,2  (-1  +  2  Vm)/(2a  -  2a  Vm)  =  -  (4.93f) 

In  summary,  to  obtain  1111,  the  desired  effective  modulus  for  this  particular  case,  the  procedure  is: 

•  Substitute  the  appropriate  numerical  values  into  Equations  4.93  and  solve  for  the 
C’s  using  an  appropriate  equation  solver. 

•  (Halculate  values  for  p. 

•  Evaluate  Till. 
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The  details  for  the  development  of  Equations  4.92  and  4.93  were  carried  out  using  Mathematica, 
and  are  to  be  found  in  Mello  (1996). 

4.7  Case  E  -  Composite  Stresses  and  Oy 

The  deformations  of  this  case  are  found  by  superposing  the  results  of  Cases  C  and  D 
examined  previously  (see  Figure  11).  Analysis  of  this  case  yields  the  effective  moduli  for 
determining  Cx  and  Oy.  This  superposition  is  analogous  to  adding  the  hydrostatic  and  deviatoric 
stresses  to  obtain  total  stress.  Recall  for  Case  C  that  the  prescribed  x  and  y  strains  are: 

e 

Ey-  2 


Cy--  2 


(4.94) 

for  Case  E.  Now  inverting  one  finds: 


e 

ex-  2 


and  for  Case  D 


e 

Ex-  2 


SO  for  the  linear  combination 


ei  62 

ex=y  +  2 


P  _  a  ^ 
2  ~  2 


61  =  Ex  +  6y 


62  =  Ex  -  Ey 


(4.95) 


The  X  and  y  composite  stresses  are  similarly  found: 
Ox  =  ai  +  02 


Oy  =  01  -  02 


(4.96) 
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Case  C  -  Change  in  Volume  Case  D  -  Cross  Section  Distortion 

with  Slip  and  Dilation 


Figure  11.-  Case  E  Deformation 


Now  from  Case  C,  Equation  4.61,  we  have 

5 

oi  =Tii  ei  +'n2e2  +  'n3ef  +  'n4  — 

M) 

and  from  Case  D,  Equation  4.89: 

<^2  =  niie2 

Substituting  Equation  4.95  into  the  above  and  putting  the  results  into  Equation  4.96,  one  obtains 
expressions  for  the  x  and  y  stresses  which  are: 

tJx  =  (Til  +  Tlll)ex  +  (Til  -  Tlii)ey  +  T12  Ez  +  TI3  Ef  +  Tl4  (4.97) 

M) 

and 

g 

Gy  =  (Til  -  Tlii)ex  +  (Tjl  +  Tlll)Ey  +  Tl2  Ez  +  TI3  Ef  +  TI4  —  (4.98) 
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which  yields  the  desired  composite  elastic  moduli.  The  remaining  effective  stress-strain  equations, 
arising  from  the  superposition  of  Cases  C  and  D,  are  found  by  substituting  Equation  4.95  into 
Equations  4.62  to  4.66  of  Case  C.  These  are: 


^z=  [Tl2ex  +  ti2ey  +  Ti5ez  +  Ti6ef  +  Ti7^  ] 


(1-p) 


1  ,  8, 

Of=  -  [Tl3ex  +  Tl3ey  +  Tl6ez  +  Tl8ef+Tl9  -] 


j  s  j  5 

T  =  d„  -  +d,2  - 


*  a  r  5 1 

C*=-2  [T\4  Ex  +T14  ey  +T17  Ez  -(-  TI9  Ef -H  TUq  -] 


,  s  ,  8 

+  <i2l  -+d22- 


(4.99) 


(4.100) 


(4.101) 


(4.102) 


This  leaves  only  one  more  composite  elastic  relation  to  be  developed. 

4.8  Case  F  -  XY  Shear 

The  final  deformation  case  considered  is  shear  in  the  xy  plane,  shown  in  Figure  12.  The 
prescribed  deformation  is  Yxy-  Recall  Case  D  which  examined  the  distortion  of  the  cross  section 
(see  Figure  lO).  Because  of  the  axisymmetric  nature  of  the  RVE,  we  may  find  the  desired  relations 
for  the  present  case  by  rotating  (transforming)  the  results  of  Case  D  by  45  degrees  (recall  the 
simple  Mohr’s  circle  type  stress-strain  transformation).  The  transformation  relations  are: 

“Cxy  =  -  (  "^2^^  )  sin(26)  +  tx'y'  cos(20) 


and 


¥  =  -  ( )  sm(2e)  +  cos(2e) 
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Figure  12.  -  Case  F,  XY  Shear  Deformation 


where  the  primed  quantities  are  from  Case  D.  For  6  of  45  degrees  and  Xx'y'  and  equal  to 
zero  as  defined  in  Case  D,  one  obtains: 


(4.103) 


Yxy  =  ex'-ey'  (4.104) 

To  find  Ox'  and  Gy'  for  Case  D  deformation  we  use  Equations  4.97  and  4.98,  now  written  as 
primed  variables 

Ox' =  (Til  +  Tlii)ex' +  (Til  -  Tiii)ey' 

Gy'  =  (Til  -  Tiii)ex'  +  (Til  +  Tiii)ey' 

which  are  substituted  into  Equation  4.103  to  get,  after  some  algebraic  manipulation: 

Xxy  =  Till  (Ex' -Ey')  (4.105) 
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Substituting  the  result  of  Equation  4.104  into  the  above  equation  leads  to  the  final  results: 

"^xy  ~  ^11  Yxy  (4.106) 

4.9  Composite  Properties 

The  results  of  Cases  A  through  E  are  now  summarized  and  organized  into  the  desired  three- 
dimensional,  homogeneous,  constituent  relations.  The  equations  that  result  are: 

Ox  =  (Til  +  Till)  ex  +  (Til  -  Till)  ey  + 112  ez  + 113  ef  + 114  ^ 

Oy  =  (111  -  Till)  Ex  +  (111  +  llll)  By  + 112  ez  +113  ef +  fi4  ^ 

1  r  S  , 

Oz=  [Tl2ex  +  il2ey+ii5ez  +  ii6ef  +  ii7  -] 

(1-p)  Tb 

1  r  S  , 

Of  =  -  [113  Ex  + 113  Ey  + 116  ez  + 118  ef +119  -  ] 
p  Id 

“Ixy  =  1111  Yxy 
'txz  =  110Yxz 
''yz  =  'n0Yyz 

,  s  ,  5 

T  =  d,l-+dl2  - 

*  r  5 

0*=  -5-  [ll4ex  +  Tl4Ey  +  Tl7Ez  +  1l9ef  +  Tlio  — ] 

^  n) 

+  d2i5^+d22^  (4.107) 
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For  implementation  into  a  finite  element  scheme,  it  is  convenient  to  write  these  equations  in  matrix 
fonn 


{a}=[D]{e} 


(4.108) 


where 


{a}T  =  {Gx,  Gy,  (l-p)Gz,  Txy,  Xyz,  txz,  pCf,  2p't,  2pG*) 


(4.109) 


£y>  Ezj  Yxy’  Yyz’  Yxz’ 


rb’ 


(4.110) 


Tll+Tlll 

111-1111 

112 

0 

0 

0 
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0 
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0 

0 

0 
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0 

0 

0 

no 

0 

0 

0 

113 

113 

116 

0 

0 

0 

118 

0 

119 

0 

0 

0 

0 

0 

0 

0 

2pdii 

2pdi2 

114 

114 

117 

0 

0 

0  , 

119 

2pd2i 

Tllo+2pd22 

(4.111) 


The  finite  element  code  developed  around  this  composite  material  model  is  limited  to  plane  strain 
conditions  and  it  is  standard  to  consider  the  x-y  plane  in  such  an  analysis.  The  above  results  are 
thus  transformed  into  this  coordinate  system  and  specialized  to  plane  strain : 

{tJ}=[D]{e}  (4.112) 


{g}T  =  {(l-p)Gx,  Gy,  Txy,  pGf,  2px,  2pG*} 


(4.113) 


{e}T=  {ex,ey,yxy,ef, 


(4.114) 


T15  Tl2  0  T16  0  T17 

Tl2  Tll+Tlll  0  T13  0  Tl4 

0  0  Tio  0  0  0 

r\6  T13  0  ns  0  T19 

0  0  0  0  2pdii  2pdi2 

-  117  Tl4  0  Tl9  2pd21  11l0+2pd22 


(4.115) 


where  the  x  direction  is  now  the  fiber  direction  (see  Figure  13). 


Y 


X 

Fiber  Direction 

Figure  13.  -  Composite  Material  Orientation  for  Plane  Strain  Analysis 
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5.  LINEAR  COMPOSITE  FINITE  ELEMENT  DEVELOPMENT 


The  composite  finite  element  was  introduced  in  Section  3.  Recall  that  the  scope  of  its 
applicability  is  limited  to  plane  strain,  linear  isotropic  fiber/matrix,  unidirectional  reinforcement, 
and  small  deformations.  For  the  finite  element  development  in  this  section,  the  bond  response  is 
also  limited  to  linear-elastic  behavior.  Bond  nonlinearity  and  its  implementation  is  presented  in 
subsequent  sections.  The  element  employs  four  degrees  of  freedom  per  node  as  discussed 
previously  (see  Figure  4).  The  composite  material  relations  developed  in  the  previous  section  are 
coupled  with  the  four  nodal  unknowns  to  form  a  specialized  finite  element  scheme. 

The  element  is  based  on  what  has  become  a  standard:  a  bilinear  displacement-based,  four- 
node  quadralateral,  two-dimensional,  isoparametric  formulation.  The  details  for  such  an  element 
can  be  found  in  many  sources  such  as  Cook  (1981)  or  Zienkiewicz  (1989).  The  steps  used  in 
deriving  this  element  will  now  be  briefly  discussed.  Details  will  be  given  where  the  element  differs 
from  the  norm  because  of  the  added  nodal  unknowns  or  its  peculiar  constitutive  relations. 

The  finite  element  equations  for  a  linear  static  analysis  can  be  posed  in  terms  of  minimizing 
potential  energy,  equating  of  virtual  work,  or  Galerkin’s  method  of  weighted  residuals,  with  the 
same  results.  This  work  will  be  extended  later  to  the  case  of  nonlinear  inelastic  material  analysis; 
therefore,  either  Galerkin’s  method  or  virtual  work  must  be  used  since  minimizing  potential  energy 
is  limited  to  elasticity.  Virtual  work  will  be  used  in  this  development  for  the  sake  of  brevity.  It  is 
recognized  that  this  method  is  less  general  than  Galerkin's  method,  however,  Zienkiewicz  has 
shown  that  for  problems  in  solid  mechanics  it  is  exactly  the  same  as  the  weak  form  of  the 
equilibrium  equations. 

The  general  governing  equation  derived  from  virtual  work  is  written  in  matrix  form  as: 


J^{Ae)T{a}dV=  IaIAuIT  {T}dA  +  {Au}T{q}dV  +  {f}  {Au}  (5.1) 


where  {T}  are  tractions,  (q)  body  forces,  and  {f}  discrete  loads.  Nodal  loads  are  the  only  loads 
considered  presently.  The  displacements  (u),  strains  (e),  and  stresses  {o}  are  defined  in  vector 
form  by: 


{u)T=  {u,  V,  Uf,  6} 

(5.2) 

{Au)T  =  {Au,  Av,  Auf,  A5) 

(5.3) 
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Yxy’  ^1’  Q2} 

(5.4) 

—  {AEx,  A£y,  Ay^^y,  A£f,  Aq],  A(J2} 

(5.5) 

{a}T  =  {(l-p)ax,  Cy,  Txy,  pGf,  2pT,  2pa*} 

(5.6) 

The  virtual  quantities  are  denoted  by  A.  The  stresses  and  strains  are  taken  from  the  constitutive 
Equation  4.107  of  the  previous  section.  Note  that  these  definitions  differ  significantly  from  that 
seen  in  a  common  homogenous  finite  element  formulation.  The  displacement  vector  contains  fiber 
displacement  and  bond  dilation.  It  is  important  to  emphasize  and  review  what  each  of  these 
vector’s  represent.  The  displacements  are: 

u  -  average  matrix  displacement  in  the  fiber  direction  (x-direction) 

V  -  average  transverse  (to  the  fiber)  displacement  (y-direction)  for  both  fiber  and 
matrix 

Uf  -  average  fiber  displacement  (x-direction) 

5  -  average  bond  dilation 

As  mentioned  previously,  this  analysis  is  limited  to  reinforcement  in  the  x-direction  only.  The 
strains  are: 

Ex  -  average  matrix  strain  in  the  fiber  direction 
Ey  -  average  transverse  strain 
Yxy  -  average  in-plane  shear  strain 
qi  -  average  slip  strain,  sA^ 
q2  -  average  dilation  strain,  5/rb 

where  the  relative  movement  of  the  fiber  and  matrix  (slip)  is 

S  =  U  -  Uf 

The  stress  quantities  are: 

(l-p)ax  -  average  matrix  stress  in  the  fiber  direction 
Oy  -  average  transverse  stress 

Txy  -  average  in-plane  shear  stress 
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pOf  -  average  fiber  stress 
2px  -  average  bond  shear  stress 
2po*  -  average  equilibrium  stress 

Recall  that  p  represents  the  fiber  volume  fraction  and  that  the  quantity  o*  is  used  to  impose 
equilibrium  at  the  bond  interface  when  dilation  is  used  as  an  unknown  in  the  composite  micro¬ 
mechanics  analysis.  This  quantity  is  set  to  zero  in  the  virtual  work  process  thus  enforcing  the 
required  equilibrium. 

For  small  deformations  the  strain  displacement  relations  may  be  written  as; 

{e}=[A]{u)  (5.7) 

or 


Ex 


ty 

Yxy 

Ef 

qi 

q2 


9x 

0 

dy 


1 


0  0  0 
d 


I 

dx 


0  0 


0 


rb 
0  0 


0  0 
0  0 


0 


rb-i 


(5.8) 


where  [A]  is  sometimes  called  the  differential  operator  matrix. 

For  a  displacement-based  finite  element  scheme,  the  body  is  approximated  as  an  assembly 
of  discrete  finite  elements  connected  at  node  points.  The  unknowns  at  the  nodes  define  the  global 
displacement  vector  (u), 

{ulT  =  {ui  VI  ufj  6i  U2  V2  . . .  5n}  (5.9) 

Within  an  element,  m,  the  displacements  are  approximated  in  terms  of  shape  functions  and  the 
nodal  displacements: 

{u}m  =  [N]n.{u}  (5.10) 
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Thus,  for  this  element  the  extra  degrees  of  fineedom,  slip  and  dilation,  are  interpolated  exactly  as  the 
x-y  displacement  components  of  a  standard  two-dimensional  homogeneous  element.  Thus, 
although  the  displacement  variables  and  constitutive  law  are  very  different  for  this  element,  only 
simple  extensions  of  an  existing  finite  element  code  arc  required  for  its  implementation. 

The  strains  and  virtual  strains  are  expressed  in  the  usual  fashion  as: 


{elm  =  [B]m  {u} 

(5.11) 

{Ae}m  =  [B]m  {Au} 

(5.12) 

where  [B]  is  the  strain-displacement  matrix  and  is  given  later. 

The  usual  virtual  work  process  leads  to  the  general  governing  equation  written  in  terms  of 
element  contributions 


nelem  nelm 

2,  f.  z  f  [N];j{,|  (5.13) 

m=l  •'''  m=l  Jvm 


where  the  only  element  loads  are  due  to  body  forces.  For  the  current  restriction  of  linear  elastic 
material  response,  one  writes 

{a}  =  [D]{e)  (5.14) 

where  [D]  is  given  in  Equation  4.15  of  the  previous  section.  Using  this  expression  leads  to 
writing  the  element  stiffness  matrix  in  the  standard  form: 

[k]m=  f  [B]jJ[D]m[B]mdv  (5.15) 


The  form  here  is  such  that  the  element  stiffness  can  be  formed  using  an  existing  finite  element 
code,  given  the  new  stress-strain  [D]  and  strain-displacement  [B]  matrices  for  the  composite 
material  and  the  four  degrees  of  freedom  (per  node)  element. 

A  general  FORTRAN  finite  element  code,  Mish  (1992),  for  two-dimensional  heat  transfer, 
HEAT2D,  was  used  as  a  basis  for  developing  this  specialized  composite  code,  hereafter  referred  to 
as  COMP2D.  HEAT2D  was  an  attractive  starting  point  as  it  contained  parameter  statements 
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governing  the  dimensions  of  all  the  global  arrays  and  came  with  user-friendly  graphical  post¬ 
processing  routines.  In  COMP2D,  as  in  HEAT2D  and  most  finite  element  codes,  the  element 
stiffness  matrix  is  calculated  in  terms  of  four  sub-matrices.  The  four  degrees  of  freedom  per  node 
of  the  composite  element  results  in  a  sub-matrix  4  by  4  in  size.  For  the  i,j  (i,j  =  1,4)  sub-matrix, 
the  elements  are  written  as 

kij=  j  [BjJlDJIBljdA  (5.16) 

The  shape  functions  are  written  in  sub-matrix  form  as: 


[N]j 


fNj  1 

Nj  0 
0  Nj 

'  NjJ 


(5.17) 


The  strain-displacement  matrix  is  found  by  multiplying  the  differential  operator  matrix  by  the  shape 
functions 


[B]j  =  [A]  [N]j 


(5.18) 


or 


[Blj  = 


aNj 

dx 

0 


dNj 

dy 

0 

Ni 

rb 

0 


0 

aNj 

dy 

aNj 

ax 

0 


0 
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0  0 
0  0 


0 

aN, 

ax 

iMi 

rb 
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0 

0 

0 

Ni 

rb  -■ 


(5.19) 
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The  stiffness  terms  for  the  i,j  sub-matrix  terms  are  now  calculated  by  doing  the  appropriate  matrix 
multiplications  in  Equation  5.16.  For  example,  one  gets 


k 


m  _ 
11  “ 


dx  dx  dy  dy 


D44+^D88]dA 


for  the  integrand  of  the  upper  left-hand  comer  term  of  this  sub-matrix.  The  integration  over  the 
element  is  performed  numerically  using  a  four-point  gauss  quadrature  scheme. 

At  this  point,  no  element  loading  due  to  either  surface  tractions  or  body  forces  is  built  into 
COMP2D,  as  the  prescription  of  discrete  nodal  loads  and  displacements  are  adequate  to  test  the 
code.  The  basic  changes  made  in  creating  CDMP2D  from  a  general  two-dimensional  heat  transfer 
code  are: 

•  Global  and  elemental  arrays  increased  to  accommodate  four  degrees  of  freedom  per 
node  as  needed. 

•  A  composite  material  algorithm  created  to  calculate  [D],  the  new  constitutive  law,  using 
the  relations  derived  from  the  micro-mechanics  analysis. 

•  A  new  16  by  16  element-stiffness  matrix  formed  using  [D]  and  the  strain-displacement 
relations,  [B],  for  the  four  degree  of  freedom  element. 

•  A  composite  stress  computation  algorithm  at  the  Gauss  points  based  on  the  new 
constitutive  law,  using  [D]. 

•  Input^output  and  format  changes. 


6.  BOND  PLASTICITY  MODEL 

The  ultimate  objective  of  this  work  is  to  implement  an  elastic -plastic  bond  model  in  a 
composite  element.  In  this  section  the  elastic-plastic  bond  model  is  described.  The  numerical 
algorithm  to  evaluate  this  analytical  model  and  the  associated  nonlinear  finite  element  scheme  are 
presented  in  Section  7. 
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A  comprehensive  plasticity  bond  model  was  developed  by  Herrmann  and  Cox  (1992),  as 
stated  previously.  The  model  can  be  applied  to  both  monotonic  and  cyclic  loading  and  is  based  on 
classic  nonassociative,  elasticity-plasticity  theory.  As  previously  noted,  the  primary  goal  of  this 
work  was  to  demonstrate  that  a  bond  model  that  includes  dilation  may  be  successfully  implemented 
in  a  composite  finite  element  scheme.  In  this  section  the  model  of  Herrmann  and  Cox  will  be 
restricted  to  monotonic  loading  and  simplified  somewhat  to  facilitate  this  initial  implementation. 
These  restrictions  may  be  removed  in  future  development.  The  complications  that  would  need  to 
be  added  impact  the  material  model  and  its  associated  numerical  algorithm  and  not  the  overall 
theory  used  in  embedding  the  bond  model  within  a  composite  representation. 

The  bond  model  has  four  main  components:  an  elastic  law,  yield  function,  hardening  law, 
and  flow  rule.  The  reader  is  referred  to  Herrmann  and  Cox  (1992)  for  the  detailed  development  of 
each  of  these  components.  The  simplified  version  used  in  this  work  is  presented  below. 

6.1  Elastic  Law 

The  bond  elastic  law  was  introduced  previously  in  Section  3  but  is  repeated  here  for 
completeness  and  is  now  written  in  terms  of  compliance: 


(6.1) 


Herrmann  and  Cox  found  the  components  of  the  compliance  matrix  to  be  proportional  to  the 
reciprocal  of  the  concrete  modulus  in  the  following  fashion: 

[A]  =  [a]/Ec  (6.2) 

where  Ec  =  concrete  modulus 

and 

an  =20.0 
ai2  =  0.86 
a2i  =  0.86 

a22  =28.6  (6.3) 
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Note  that  these  values  are  twice  those  reported  by  Herrmann  and  Cox,  due  to  the  fact  that  in  their 
work  the  slip  and  dilation  strains  were  normalized  with  respect  to  the  bar  diameter  as  compared  to 
the  bar  radius  here. 

6.2  Yield  Function 

Consistent  with  classical  plasticity  models,  a  yield  function  or  surface  is  used  to  define  the 
limits  of  elasticity.  For  the  behavior  of  concrete  bond,  this  function  must  harden  and  soften  as 
damage  occurs.  The  yield  function  of  Herrmann  and  Cox  simplifies  to  the  form: 

F=  i -C(D)a,(^)“2  (6.4) 


where 


Ft  =  concrete  tensile  strength 

C(D)  =  hardening  law  which  is  a  function  of  damage,  D  (to  be  discussed  subsequently) 
ai,  a2  =  parameters  introduced  by  Herrmann  and  Cox 


This  function  is  plotted  in  Figure  14  for  the  case  of  no  damage,  which  should  be  compared  to 
Herrmann  and  Cox's  Figure  19.  This  power  law  yield  surface  is  a  simplification  of  the  blended 
power  law-exponential  function  used  by  Herrmann  and  Cox.  This  simplification  is  made  in  part  to 
allow  for  the  use  of  a  nonassociative  flow  rule  that  is  simply  related  to  the  yield  surface.  A 
discussion  of  this  follows  the  definition  of  the  flow  rule.  Because  of  this  simplification,  the 
fidelity  of  the  predictions  to  experimental  results  is  less  than  that  achieved  by  Herrmann  and  Cox 
when  there  is  significant  damage  combined  with  low  confinement  (bond  normal  stress). 

The  damage,  D,  serves  as  the  internal  variable  governing  the  evolution  of  the  yield  surface 
and  the  resulting  plastic  flow.  Herrmann  and  Cox  found  the  best  damage  measure  to  be: 


D  = 


q?  Tb 


Sl  cos<|) 


(6.5) 
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U/1 
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Function 


where 


Qj  =  plastic  slip  strain 

rb  =  reinforcement  radius 

Sl  =  reinforcement  lug  spacing 

(j)  =  reinforcement  lug  angle 

Thus,  damage  is  essentially  a  function  of  the  plastic  slip. 

6.3  Hardening  Law 

The  hardening  and  softening,  or  changing  size,  of  the  yield  surface  is  controlled  by  the 
value  of  the  C(D)  term  in  Equation  6.4.  The  hardening  law  of  Herrmann  and  Cox  was  given  as  a 
differential  equation  in  damage,  D.  This  is  simplified  in  this  work  to  the  following  functional 
relation: 


C(D)  =  03  +  04  (1  -  e-“5D)  -  06  (1  -  e-“7D)  (6.6) 

This  is  accomplished  by  numerically  solving  the  ordinary  differential  equation  of  Herrmann  and 
Cox  and  fitting  the  parameters  03  -  07  to  best  match  their  results.  This  was  accomplished  using 
Mathematica.  The  details  of  this  process  can  be  found  in  Mello  (1996).  Figure  15  compares  the 
two  results. 

6.4  Nonassociative  Flow  Rule 

The  nonassociative  flow  rule  which  governs  the  relation  between  plastic  slip  and  plastic 
dilation  strains  was  defined  by  Herrmann  and  Cox  in  terms  of  function  g,  where: 


g  = 


(6.7) 
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Simplified  Function  Law 


H 
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Figure  15.  -  Bond  Model  Hardening  Law 


Herrmann  and  Cox  found  g  to  be  a  complex  function  of  bond  normal  stress,  damage,  plastic 
dilation,  and  nearly  a  dozen  parameters.  Once  again  this  relation  is  simplified.  A  less  complex 
flow  rule  is  derived  by  rotating  back  a  prescribed  constant  angle  0r  from  the  normal  to  the  yield 
surface  as  defined  by  Equation  6.4  (the  angle  is  a  calibrated  quantity).  This  is  shown  in  Figure  16. 
Using  this  assumption 


•pi 

qf 


=  tanG 


(6.8) 


where  0  is  the  angle  between  the  flow  direction  and  the  x  axis 

0  =  0*  -  0r 


and 


^  ,  3F  .  ap  , 

0*  =  tan-i  ( —  /  —  ) 

do  dx 


letting 


.  ap  ,  ap 

do  dx 

it  is  easy  to  show  that: 


(6.9) 


(6.10) 


(6.11) 


h  -  tanGr 
1  +  h  tanGr 


(6.12) 


The  only  parameter  to  be  determined  in  this  case  is  the  rotation  angle,  0r,  assuming  the  yield 
surface  and  its  associated  derivatives  have  been  defined. 

The  calibration  of  0r  is  done  using  the  simplified  yield  surface  Equation  6.4,  its  derivatives, 
the  damage  function  C(D)  of  Equation  6.6,  and  the  experimental  results  of  Gambarova  as  reported 
in  Figure  30  of  Herrmann  and  Cox  (1992).  The  angle  was  adjusted  so 
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that  the  results  for  g  best  match  those  of  Herrmann  and  Cox  given  in  their  Figure  25.  Details  are 
given  in  Mello  (1996).  The  best  fit  is  found  to  be: 

0r  =  0.21  rad  (6.13) 

The  analytical  yield  function,  hardening  law,  and  flow  rule  are  implemented  in  a  numerical 
nonlinear  material  algorithm  that  evaluates  the  bond  stress  state,  given  a  strain  history.  Details  of 
this  algorithm  and  its  implementation  in  the  global  composite  finite  element  scheme  follow  in  the 
next  section. 


7.  NONLINEAR  BOND  FINITE  ELEMENT  IMPLEMENTATION 

The  incorporation  of  the  simplified  elastic-plastic  bond  model  into  the  composite  finite 
element  scheme  is  detailed  in  this  section.  First,  the  finite  element  solution  of  nonlinear  plasticity 
problems  is  discussed  in  general  terms.  The  process  is  outlined  briefly  to  understand  what  is 
required  of  the  bond  model  and  composite  material  algorithm.  There  are  many  references  in  the 
area  of  nonlinear  finite  element  analysis.  Those  used  most  here  are  Zienkiewicz  and  Taylor 
(1991),  Owen  and  Hinton  (1980),  and  Bathe  (1982). 

7,1  Nonlinear  Finite  Element  Solution 

The  source  of  nonlinearity  for  this  study  is  limited  to  the  elastic-plastic  behavior  of  the 
bond.  It  is  recognized  that  the  behavior  of  the  matrix  and  reinforcing  fibers  may  be  nonlinear  also, 
but  in  this  study  the  goal  is  the  initial  implementation  of  a  plasticity  bond  model  including  dilation 
effects.  Because  the  analysis  also  assumes  small  deformations,  no  reformulation  of  the  strain- 
displacement  relations  is  required. 

An  incremental  Newton-Raphson  iterative  method  is  used  to  solve  the  nonlinear  problem. 
As  in  the  linear  case,  the  problem  is  formulated  in  terms  of  unknown  nodal  displacements.  The 
global  residual  is  defined  for  a  load  increment  as: 

{R(u„+i))  =  {P(un+l)}  -  {fn+ll  (7.1) 

where  {f}  are  the  applied  forces  and 
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{P}=J^[B]T{a}n+i  dv 


(7.2) 


is  the  vector  of  internal  forces.  The  vector  {P}  is  a  nonlinear  function  of  displacement  as  a  result 
of  the  bond  nonlinearity.  An  iterative  solution  to  the  problem  is  sought  and  found  when  the 
residual,  (R),  of  Equation  7.1  is  near  zero.  For  a  Newton-Raphson  iterative  scheme.  Equation 
7.1  is  written  as  a  truncated  Taylor  series  and  set  to  zero,  making  the  residual  small: 


(R^i))  -  {R(u„' ,))  f  I  ^  1  sui:;  =0 


(7.3) 


where 


r  — 1  ‘ 

^  9u  ^n+l 


is  the  global  Jacobian.  Note  the  superscripts  are  the  iteration  counters.  The  iterative  correction  is 
desired  and  is  found  by  solving 


m„i.(8uii:i=-{R)  ‘ 


(7.4) 


which  is  identical  in  form  to  the  system  of  equations  that  occurred  for  the  linear-elastic  case.  The 
displacement  approximation  is  updated  each  iteration  in  the  following  manner. 


(u) 


i+1 

n+1 


{u}„;i  +  {6u) 


i+1 

n+1 


(7.5) 


where  Un  is  the  displacement  at  the  end  of  the  last  load  increment.  The  Newton-Raphson  iterative 
process  is  repeated  until  convergence,  that  is,  until  some  norm  of  the  residual  is  sufficiently  small. 
With  this  procedure,  most  of  the  structure  of  a  linear-elastic  finite  element  code  is  used.  The 
important  new  components  on  the  global  level  are  the  global  residual  and  Jacobian. 
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7.1.1  Global  Residual.  Recall  the  global  residual  is: 


{R}=J^^j[B]T{a}dv  -  {f) 


This  requires  knowing  the  composite  stress  state  {a}  for  the  given  displacement  estimate  and 
resulting  strains.  Obtaining  the  stresses  for  this  composite  scheme  is  a  two-step  process: 

1.  A  suitable  numerical  algorithm  is  used  to  determine  bond  model 
stresses  for  a  given  bond  strain  history. 

2.  The  material  micro-mechanics  model  is  then  used  to  determine  the 
composite  stresses  given  the  bond  stresses  and  the  composite  strains. 


The  second  step  is  addressed  first,  as  it  is  simple  compared  to  the  task  of  determining  the  bond 
stresses  in  Step  1. 

The  composite  stress  relations  for  a  nonlinear  bond  stress  model  arc  obtained  by  returning 
to  deformation  Case  C  of  Section  4  which  involves  the  bond.  Instead  of  substituting  the  linear- 
elastic  bond  stiffness  in  the  constituent  virtual  work  expression  for  this  deformation  mode,  it  is  left 
in  terms  of  the  bond  stresses.  Equation  4.66  is  now: 


o*=  ^  [Tl4ei+Ti7e2  +  Ti8ef  +  Tlio^]+o 


(7.6) 


The  virtual  work  of  the  composite  and  constituents  is  equated  as  before  but  now  the  composite 
stress-strain  relations  that  result  are  in  terms  of  the  bond  stresses  as  well  as  the  strains.  These 
resulting  composite  stress-strain  relations  for  a  nonlinear  bond  stress  model  are: 


1  r  S  , 

Ox=  [Tl5ex+Tl2ey+Tl6ef  +  Tl7  :^] 


(7.7) 


=  712  Ex  +  (T|l  +  Till)  Ey  +  Tl3  Ef  +  T14 


rb 


(7.8) 


'txy  -  TlO  Yxy 


(7.9) 
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(7.10) 


1  8 
0{=  -  [Tleex  +  TlcEy  +  TlgEf  +  Tla  “1 


T=t 


(7.11) 


a*  = 


^  [mex+Tl4ey  +  Ti9ef  +  Tiio -]  +  o 


(7.12) 


where  x  and  o  are  the  nonlinear  bond  stresses  determined  in  Step  1  from  a  numerical  evaluation 
algorithm  for  the  bond  model. 

The  first  step  in  determining  the  bond  stresses  turns  out  to  be  quite  difficult  due  to  the 
nonassociative  plasticity  and  a  yield  function  that  is  not  defined  for  tensile  stresses.  Before 
developing  an  algorithm  to  determine  the  bond  stresses  from  the  bond  model,  we  must  determine 
what  additional  information  the  global  solution  scheme  requires  from  such  an  algorithm. 


7.1.2  Global  Jacobian.  Recall  the  global  Jacobian  equation  was  defined  as: 


[J] 


which,  for  constant  loads  and  small  deformations  is  written  as: 

J=f  [B]T[^][B]dv  (7.13) 

Jv  de 

The  global  Jacobian,  [J],  which  can  be  thought  of  as  a  consistent  stiffness  matrix  in  the  tangent 
direction,  is  built  up  from  the  global  material  Jacobian  or  consistent  tangent  material  stiffness  and 
the  linear  strain-displacement  matrix.  The  global  material  Jacobian  for  the  composite  is  found  by 
substituting  the  bond  material  Jacobian 


Jbg  = 


9x  3t 
8qi  9q2 


L  9qi  3q2  J 


(7.14) 
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into  the  composite  material  relations  in  place  of  the  elastic  bond  matrix.  Thus,  the  two  items 
required  and  yet  to  be  derived  are  the  bond  stresses  and  the  bond  material  Jacobian. 

7.2  Numerical  Implementation  of  the  Bond  Model 

To  complete  the  implementation  of  the  nonlinear  bond  model,  a  numerical  material 
algorithm  is  required.  This  algorithm  must  use  the  yield  function,  hardening  law,  and  flow  rule  for 
the  bond  model,  as  given  in  Section  6,  to  calculate  the  bond  stress  state  and  material  Jacobian  for  a 
given  strain  history. 

7.2.1  Cutting'Plane  Algorithm.  The  cutting-plane  algorithm  of  Hughes  (1987)  was 
initially  selected  and  used.  This  procedure  uses  a  return  mapping  algorithm  to  integrate  the 
inelastic  constitutive  equations.  The  algorithm  worked  adequately  for  problems  where  significant 
bond  confinement  pressure  existed  but  had  numerical  difficulties  with  more  general  bond  states. 
The  trouble  resulted  from  the  fact  that  the  algorithm's  iterative  process  employed  an  elastic  trial 
stress.  Using  the  bond  elastic  law,  the  trial  stress  (when  no  initial  confinement  existed)  resulted  in 
a  positive  bond  normal  stress.  Positive  normal  stresses  are  not  defined  mathematically  or 
physically  for  the  bond  model.  A  variety  of  remedies  were  attempted.  The  most  promising  used  a 
new  yield  surface  with  a  fictitious  extension  of  the  surface  in  the  positive  stress  region  blended  to 
the  desired  power  law.  This  is  shown  in  Figure  17.  This  fix  worked  for  most  situations  but  had 
stability  problems  in  the  case  of  little  or  no  confinement,  which  occurs  at  the  start  of  most  practical 
problems. 

The  solution  to  the  above  dilemma  was  the  development  of  a  new  more  robust  algorithm. 
The  algorithm  is  based  on  recent  work  by  Herrmann  (1996).  His  Reduced-Newton  algorithm  had 
been  used  successfully  on  complicated  soil  material  models  that  were  mathematically  similar  to  the 
bond  model.  The  details  of  the  application  of  this  algorithm  to  the  bond  model  follow. 

7.2.2  Reduced-Newton  Material.  The  basic  idea  of  this  algorithm  is  to  analytically 
reduce  the  model  equations  to  a  minimum  number  of  equations  and  unknowns.  In  this  case,  two 
pair  of  nonlinear  ordinary  differential  equations  were  obtained.  These  equations  are  expressed  in 
terms  of  the  bond  stresses  which  are,  in  turn,  functions  of  a  pseudo-time  variable,  t,  representing 
the  history  of  the  strain  increment.  These  equations  are  solved  using  a  fully  implicit  (or  backward 
difference)  time-marching  scheme,  and  because  the  equations  are  nonlinear,  a  Newton-Raphson 
iteration  as  well. 
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Figure  17.  -  Yield  Function  Including  Fictitious  Bond  Tension 


Material  Algorithm  Inout/Output 


Input  to  the  algorithm  is  derived  from  the  solution  for  the  previous  solution  step  and  the 
current  estimate  of  the  displacement  increment: 


{Aqi  Aq2} 

-  estimate  of  strain  increment  for  current  step  (seeking  stress  at  the 
end  of  this  increment) 

{  qib  q2b) 

-  total  strain  at  the  end  of  the  last  step 

fqi?  1 

-  plastic  slip  strain  at  the  end  of  the  last  step 

{'Cb  Ob) 

-  total  stress  at  the  end  of  the  last  step 

Required  output  from  the  material  algorithm  is: 

{'^n  <yn} 

-  new  stress  state  (to  be  used  in  the  global  residual) 

(11?  1 

-  new  plastic  slip  strain 

UbgI 

-  bond  material  Jacobian 

Time  Stepping 

The  application  of  the  given  strain  increment  is  defined  in  terms  of  a  pseudo-time  variable 
within  the  increment  or  step  as  shown  in  Figure  18.  The  quantities  known  at  the  beginning  of  the 
increment,  or  t=0,  are  designated  by  the  b  subscript  while  the  n  denotes  those  to  be  found  at  the 
end  of  the  increment  (t=l).  Thus,  the  strain  over  the  given  increment  is  written  as  a  function  of 

time, 


qi(t)  =  qib  +  Aqi  t 

q2(t)  =  q2b  +  Aq2  t  (7.15) 
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t,  pseudo  time  variable 

1.0 


Finite  Difference  Sub-Step  Scheme 


Figure  18.-  Strain  Increment  Time  Variable  and  Sub-Stepping 


with  the  underlying  assumption  that  the  strain  increments  are  sufficiently  small  so  that  the  strain 
histories  are  proportional  over  the  increment.  Now  one  gets  the  strain  rates  by  taking  the  time 
derivatives  of  the  above  expressions. 


dqi  _ 
dt  ■ 


Aqi 


dq2 

dt 


=  Aq2 


(7.16) 


The  dot  (•)  notation  is  employed  so  that  the  rates  are  written: 

qi=Aqi  q2  =Aq2  (7.17) 

The  development  that  follows  is  limited  to  plastic  loading,  that  is  the  initial  stress  state,  designated 
by  b,  is  initially  on  the  yield  surface  and  a  strain  increment  results  in  plastic  flow,  evolution  of  the 
yield  surface,  and  resulting  stresses  on  a  new  yield  surface. 

Elastic  Law 

The  bond  model  elastic  law  given  in  Equation  6.1  is  written  in  equation  form  as: 
q®  =AiiT  +  Ai2a 

(7.18) 

q2  =  Ai2  'T  +  A22  o 

Plastic  Strain 

The  elastic  strain  rates  are  found  by  differentiating  the  above  expression  to  get 
q®  =AiiX  +  A120 

(7.19) 

q2  =  Ai2  X  +  A22  o 
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During  yielding,  the  material  will  behave  partly  elastic  and  partly  plastic.  Thus,  during  any 
increment  of  stress,  the  changes  in  strain  are  assumed  to  be  divisible  into  elastic  and  plastic  parts, 


•  *6  ,  ‘pl 

qi=qi+  qi 


•  *6  .  ‘pl 

q2=q2+  q2 


(7.20) 


One  then  solves  for  the  plastic  strain  rates  by  substituting  the  elastic  rates  which  are  given  in 
Equation  7.17, 


qP*  =Aqi  -  All  X  -  Ai2a 

(7.21) 

q^*  =  Aq2  -  Ai2  X  -  A22  O 


Hardening  Law 

Since  the  hardening  law  that  controls  evolution  of  the  yield  surface  is  in  terms  of  the  total 
plastic  slip,  it  is  necessary  to  integrate  the  above  rate  equations. 


I 

I 


i 


I 


I. 


■  <  d,-=  q?  (,)-,,?■ 


(7.22) 


Performing  this  integration  and  solving  for  the  plastic  slip  at  time,  t,  results  in  the  following: 


qf  (t)  =  qi?  +  Aqit  -  Aii(t(t)  -  xb)  -  A12  (a(t)  -  Ob)  (7.23) 

The  hardening  law  is  finally  written  as: 


CD(t)  =  a3  +  04  (1 


.  e““5^  q?  ) 


(7.24) 
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where 


Dn  = 


rb 


(Slcos<1)) 


(7.25) 


is  the  coefficient  that  normalizes  plastic  slip  to  obtain  damage. 
Yield  Function 


The  yield  function  of  Equation  6.4  is  rearranged  algebraically  and  is  treated  as  a  residual  Ri 
in  the  numerical  algorithm, 


Ri  =  (#)'^  +  (CD(t)ai)P  (^) 


(7.26) 


The  p  power  is: 


P=- 

012 


(7.27) 


This  form  of  the  yield  surface  is  numerically  the  same  as  given  in  Equation  6.4  but  has  the 
advantage  that  it  may  be  evaluated  for  initial  values  of  t  and  o  that  are  zero,  which  is  the  case  for 

most  real-world  problems.  This  equation  is  cast  as  a  residual,  thus,  if  the  unknowns,  stress  at  time 
t,  are  such  that  they  converge  to  the  yield  surface,  then  this  residual  is  zero. 

Flow  Rule 

The  flow  rule  of  Equations  6.7,  6.11,  and  6.12  is  treated  as  a  second  residual  and  is 
written  as: 


h(t)+tan6r 
A-h(t)tan0r  qP‘(t) 


(7.28) 


where 
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h(t)  = - ^ ^  (7.29) 

(CD(t)ai)P 

This  form  of  the  residual  had  division  by  zero  problems  in  some  cases  and  thus  was  rearranged  so 
as  to  be  mathematically  well  defined: 

R2  =  92  (h(t)  +  tanGr)  -  qf  (1  -  h(t)  tan0r)  (7.30) 

The  plasticity  equations  of  Section  6  thus  reduce  to  Equations  7.26  and  7.30  with  the  appropriate 
substitutions.  These  form  a  coupled  set  of  nonlinear,  ordinary  differential  equations  in  the 
unknown  stresses  T  and  a  that  are  solved  as  functions  of  the  pseudo-time  variable  t. 

Finite  Difference  Scheme 


A  finite  difference  method  is  used  to  solve  the  coupled  nonlinear,  ordinary  differential 
equations  that  define  the  plasticity  model.  To  accomplish  this,  the  time  step  is  divided  into  uniform 
sub-steps,  as  shown  in  Figure  18.  Using  a  backward  difference  scheme  (Kreyszig,  1983),  the 
unknown  stresses  and  associated  first  derivatives  are  written  as: 

t(t)  =  Tn  o(t)  =  an 

X  =  - 

Atn 


Gn-On-l 

G=  - 

Atn 


(7.31) 


These  difference  relations  are  substituted  into  the  equations  that  make  up  the  residuals  Ri  and  R2, 
the  coupled  ordinary  differential  equations  of  yield  and  flow.  The  result  is  a  pair  of  nonlinear 
algebraic  equations  for  Xq  and  On. 
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Newton-Raphson  Iteration 


A  Newton-Raphson  iterative  scheme  is  used  to  solve  the  nonlinear  algebraic  equations  for 
each  time  step.  The  unknown  stresses  for  this  process  are  written  in  iterative  fonn  for  time,  tn,  as: 


(7.32) 


Expanding  the  residuals  in  a  truncated  Taylor  series  and  making  the  (I)  residuals  small,  one  gets  in 
the  normal  Newton-Raphson  fashion. 


^  =  *^1(1-1)  ^  ^"(I)  ^"(1) 

^  ^  (I-l)  ^  ’  dOn  (i-i)  ® 


(7.33) 


0  =  ^"(1)+^  ^no) 

5Xn(I-l)  3an(l.i) 


which  results  in  a  system  of  equations  for  the  Newton-Raphson  stress  corrections: 


~  aRi  aRi  ~ 

dXn  30n 


8R2  3R2 


(I-l) 


8Xn 

80n 


"1 
R2  Ja-1) 


(7.34) 


or 
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(7.35) 


■  Jll  Jl2  ■ 

-  hi  hi  - 


8xn 

SOn 


The  Jij  is  the  Jacobian  used  for  solving  for  the  stress  corrections  above. 

The  first  residual  equation  is  now: 

Ri  =  (^)P  +  (CD„ai)l*(|^)  (7.36) 

where 

Cdh  a3  +  04  (1  -  e““5DN  ‘ifi )  -  06  (1  -  e-“7DN  qfi )  (7.37) 

=qi6*  +  Aqitn- Aii(Xn-'tb)-Ai2(ctn-<5b)  (7.38) 

The  Ill  and  J12  Jacobian  terms  are  found  by  taking  the  derivatives  of  this  first  residual  equation. 
For  the  Jll  term. 


Jll  = 


aRi 

axn 


where 


|-  ( ^  )(P-1)  +  pai(CD„  ai)(P-l) 


dxn 


^  =  [a4a5e"“5DN‘‘i?  -  abave""^!^!?  ]  ^ 

aXn 


(7.39) 


(7.40) 


(7.41) 
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(7.42) 


pi 


9Xn 


=  -Aii 


For  the  J12  term. 


where 


where 


Ji2=  r* 

(7.43) 

(7.44) 

"  =  [a4a5e  “5*^N^n  -  oeaTe  ]  - 0- 

aOn 

(7.45) 

A 

...  --A12 

aon 

(7.46) 

The  second  residual  equation  is: 

^2  =  T2n  ‘  tanGf) 

(7.47) 

hn  =  |3(^)*‘'’/(CD„ai)l’ 

(7.48) 

Atn  At„ 

(7.49) 

qSl=A<12-Al2(^)-A22(5E5tl) 

Atn  Atn 

(7.50) 

The  associated  Jacobian  terms  J21  and  J22  are  now  found.  For  the  J21  term, 
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(7.51) 


8R2 

dt 


where 


3R2  _  9hn 
3Xn  dXn 


(hn  +  tanOr) 


-  -r^  (1  -  hn  taner) 
atn 


(7.52) 


^hn 

dtn 


(Cd„  ai#  ^  )<P-'>  oi  (CD,  ai#-»  ^ 

(  (CDn  ai)P)^ 


(7.53) 


^qfn  _  All 
dXn  Atn 

_  Ai2 
dXn  Atn 

For  the  J22  term, 


h2  = 


aR2 

aon 


(7.54) 


(7.55) 


(7.56) 
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where 


9R2  _  3hn 


aqP' 


( tan0r  +  8q^J[ )  +  — ^  (hn  +  tanOf) 


aon 


dqP* 


In 


8an 


(1  -  hn  +  tanGr) 


^  =  -  p2  ( 
8<Tn 


a.  (Cd„  ^ 

8an 


^fn  _  An 
30n  Atn 

^q2n  _  A22 

8cJn  Atn 


(7.57) 


(7.58) 


(7.59) 


(7.60) 


These  equations  are  used  to  solve  for  tn  and  On  using  Newton-Raphson  iteration  for  each  time 
step. 

Bond  Material  Jacobian 


Recall  that  in  addition  to  the  new  bond  stresses,  the  global  finite  element  routine  requires 
the  bond  material  Jacobian  in  order  to  form  the  global  Jacobian.  The  material  Jacobian  is: 


[Jbg]  = 


1 _ 

8qi  8q2 

r  Jgii  Jgi2  ■ 

do  8o 

L  JG21  JG22  - 

-  5qi  8q2  - 

(7.61) 


which  for  the  particular  increment  of  strain  is  written: 
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dXn  9Xn 


dAqi  dAq2 
-  9Aqi  9Aq2  - 


(7.62) 


To  derive  these  Jacobian  terms  we  return  to  the  two  residual  equations  of  the  yield  surface  and  the 
flow  rule,  Equations  7.36, 7.47.  Now,  for  the  assumed  plastic  loading, 

aAqi 


^=0 

9Aq2 


(7.63) 


that  is,  the  stress  state  remains  on  the  yield  surface  for  an  increment  of  strain.  Likewise, 


3R2 

Aqi 


=  0  , 


iS2-=o 

9Aq2 


(7.64) 


the  flow  rule  must  be  satisfied. 

Implicit  partial  differentiation  of  the  residual  equations  and  the  chain  rule  are  used  to  form 
Equations  7.63  and  7.64: 

^  Rln  /  ■.  ^  9Rln  /  \  _  q 

9Aqi  dtn  9Aqi  dOn  9Aqi 

^Rin  I  Rln  (-^)+  ^  (i5n_)  =  0 
9Aq2  dXn  9Aq2  dCn  9Aq2 

9R2n  ^  ^  /  ^Xn  ,  9R2  .  ._q 

dAqi  0Tn  9Aqi  9an  9Aqi 

^+E2a(il!L)^,^(i2!L)  =  o  (7.65) 

3Aq2  9Xn  9Aq2  90n  9Aq2 
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The  results  are  four  simultaneous  equations  in  the  desired  global  material  Jacobian  terms.  Note 
that  some  of  these  terms  are  identical  to  the  Jacobian  of  Equation  7.35.  Additionally,  the 
derivatives  of  the  residual  equations  with  respect  to  the  strain  increments  Aqi  and  Aq2  need  to  be 
derived.  The  four  new  terms  that  result  are: 


aRi 

aAqi 


=  Oil  P(CDn  ai) 


(p-1)  ^CDn 

aAqi 


^  =  [0405  D„e-“5l>N<ll? .  0(6a7D„e-“7‘>N<'l5'  ]  ^ 

aAqi  aAqi 


aAqi 


and 


=0 

aAq2 


(7.66) 


(7.67) 


(7.68) 


(7.69) 


and 


^  ^  (qi;taner+q|;)+  (h„  +  lane,) 


aAqi  aAqi 


aAq2 


(7.70) 


aAqi 


(1  -  hn  tanOr) 


aAqi  *"1  aAqi 


j  ^  A^  a^n-i  ^  Ai2  aOn-i 
aAqi  At  aAqi  At  aAqi 


(7.71) 


(7.72) 
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pi 

°^2n  ^  j  ^  ^  ^tn-1  ^  A22  ^gn-1 
9Aq2  At  0Aqi  At  9Aqi 


(7.73) 


Note  that  the  derivatives  of  the  stresses  at  the  previous  time  step  will  be  known  for  the  present  time 
step  for  which  these  equations  are  formulated.  The  final  term  is: 


aR2  ^ 

^^2  8Aq2 


(hn  +  tanGr) 


(l-hntan0r) 

dAq2 


(7.74) 


^  ^Vl  ^  ^  ^^n-l 
3Aq2  At  Aq2  At  dAq2 


(7.75) 


pi 

^^2n  _  j  ^  Ai2  9tn-l  ^  A22  ^^n-l 
9Aq2  At  3Aq2  At  3Aq2 


(7.76) 


All  the  terms  required  to  calculate  the  bond  material  Jacobian  have  been  found.  The  system  of 
equations  for  the  Jacobian  is  rewritten  as: 

ail)JG|l+<]l2)]G21=  ^ 

dAqi 

(J21)JGii  +  (J22)JG21=^  (7-77) 

oAqi 

(Jll)  JGi2  +  (J12)  JG22  =  0 

(J21)  JGi2  +  (J22)  JG22  =  ^  (”7.78) 

aAq2 

The  resulting  two  pair  of  equations,  each  in  two  unknowns,  are  solved  numerically  using  Cramer’s 
rule.  Note  the  whole  process  of  calculating  the  bond  material  Jacobian  is  not  performed  until  the 
converged  stresses  are  determined  by  appropriate  sub-stepping.  This  is  accomplished  by  saving 
the  converged  stress  history.  Later,  the  recursive  bond  material  Jacobian  equations  are  solved 
without  iteration  given  the  converged  stress  history.  This  saves  substantial  computation,  as  this 
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process  must  be  repeated  for  each  global  Newton-Raphson  iteration  and  for  each  quadrature  point 
in  the  finite  element  model.  Details  of  the  FORTRAN  code  used  to  evaluate  the  model  can  be 
found  in  Mello  (1996). 


8.  ELEMENT  VERIFICATION  TESTS 

The  element,  its  material  model,  and  FORTRAN  code  implementation  were  tested  during 
various  stages  of  development.  The  three  main  test  sets  were  linear  composite  element  tests,  bond 
model  tests,  and  structural  tests.  Examples  of  each  are  discussed  below. 

8.1  Linear  Composite  Element  Tests 

This  stage  of  testing  was  performed  after  development  of  the  composite  material  model  and 
the  creation  of  the  four  degree  of  freedom  linear-elastic  element  based  on  the  composite  material 
relations.  The  analyses  were  aimed  at  testing  the  micro-mechanics  material  model  accuracy  and 
verifying  that  the  finite  element  stiffness  matrix  and  global  arrays  were  modified  to  handle  four 
degrees  of  freedom  correctly.  The  results  reported  here  are  for  no-slip  micro-mechanics,  uniform 
slip,  and  shear-lag  problems. 

8.1.1  No-SIip  Tests.  The  no-slip  micro-mechanics  tests  were  performed  to 
demonstrate  that  the  composite  material  model  would  specialize  to  the  case  of  perfect  bonding 
between  fiber  and  matrix,  that  is  the  case  of  no-slip.  This  was  done  in  the  four  degrees  of  freedom 
finite  element  setting  using  a  penalty  formulation  (Cook,  1981)  to  constrain  the  fiber  displacement 
to  that  of  the  matrix  (no-slip).  Additionally,  the  dilation  was  set  to  zero  for  all  nodes  in  the  model. 
Single  and  multiple  element  tests  were  performed.  Three  of  the  tests  were  uniform  displacement  in 
the  fiber  direction,  the  transverse  direction,  and  a  state  of  shear  in  the  xy  plane  (plane  of  the  fibers). 
These  tests  are  essentially  analytical  versions  of  the  three  basic  orthotropic  composite  material 
characterization  cases.  The  results  were  compared  to  a  perfect  bonding  micro-mechanics  material 
model  based  on  a  two-dimensional  finite  element  analysis  of  a  hexagonal  packed  fiberglass/epoxy 
composite  (Schulz,  1965).  A  comparison  of  results  is  presented  below  in  Table  1.  The  small 
differences  are  due  to  the  approximation  of  modeling  the  hexagonal  cross  section  of  the  unit  cell  as 
circular  (see  Figure  5). 
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Table  1.  No-Slip  Micro-Mechanics  Test  Results 


Strain  Fiber  Direction  (x) 

Method 

Transverse  Strain 
(in.) 

Fiber  Direction  Stress 
(psi) 

Hex  pack,  micro-mechanics 

-0.00444 

87660 

Composite  finite  element 

-0.00457 

88030 

Percentage  difference  (%) 

2.9 

0.42 

Strain  Transverse  Direction  (y) 

Method 

Fiber  Dir.  Strain  (in.) 

Transverse  Stress  (psi) 

Hex  pack,  micro-mechanics 

-0.00145 

28630 

Composite  finite  element 

-0.00146 

28080 

Percentage  difference  (%) 

0.70 

1.9 

In-Plane  Shear  (xy)  Fiber-Transverse  | 

Method 

Composite  Shear  Stress  (psi) 

Hex  pack,  micro-mechanics 

7150 

Composite  finite  element 

7043 

Percentage  difference  (%) 

1.5 

. - . -- .  - 

These  results  show  the  correctness  of  the  composite  micro-mechanics  model.  They  also 
show  that  the  finite  element  and  its  embedded  linear  micro-mechanics  material  model  can 
specialize  to  the  case  of  perfect  bonding. 

8.1.2  Uniform  Slip.  The  uniform  slip  case  tested  the  implementation  of  the  bond 
elasticity  model  by  isolating  its  response.  These  analyses  included  one-  and  four-element  models 
where  the  transverse,  Uy^  and  matrix,  Ux,  displacements  were  set  to  zero  while  a  fiber 

displacement  or  slip  was  prescribed.  For  prescribed  slip,  the  dilation  and  bond  stresses  that  result 
were  compared  to  simple  analytical  calculations  based  on  the  bond  elasticity  relations.  A 
comparison  of  results  is  given  in  Table  2.  The  linear  bond  model  is  implemented  correctly  as 
shown  by  these  results. 
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Table  2.  Uniform  Slip  Test  Results 


Method 

Dilation 

Bond  Shear  Stress 

Bond  Normal  Stress 

Analytical 

5.290x  10'^ 

6652 

-6649 

Finite  Element 

5.291  X  lO'^ 

6653 

-6649 

8.1.3  Shear-Lag  Tests.  The  shear-lag  test  case  compared  a  classical  shear-lag  solution 
to  finite  element  results.  The  test  and  corresponding  composite  finite  element  model  is  shown  in 
Figure  19.  This  test  simulated  pulling  reinforcement  from  a  fully-fixed  prism  of  concrete. 

An  approximate  analytical  solution  was  found  by  reducing  the  problem  to  an  ordinary 
differential  equation  in  slip.  Figure  20  presents  a  comparison  of  the  analytical  and  finite  element 
results.  The  composite  element  is  able  to  reproduce  linear  shear-lag  behavior  that  is  associated 
with  an  edge  effect  with  a  very  simple  mesh. 

8.2  Bond  Model  and  Algorithm  Tests 

Tests  were  performed  on  the  bond  model  and  its  associated  numerical  algorithm  prior  to 
implementation  into  the  composite  finite  element.  Testing  was  performed  by  driving  the  bond 
model  material  algorithm  with  what  amounts  to  a  one-gauss  point  finite  element  scheme  written 
in  terms  of  strains  and  stresses.  This  driver  was  patterned  after  a  similar  one  developed  by 
Herrmann  and  Kaliakin  (1987)  that  was  used  to  evaluate  soil  material  models  and  algorithms. 

Results  of  two  test  cases  are  presented.  These  results  correspond  to  some  of  the  bond 
material  characterization  tests  that  Herrmann  and  Cox  (1992)  used  in  development  of  their  bond 
model.  The  goal  of  the  present  tests  is  to  make  sure  that  the  simplified  bond  model  gives  results 
that  qualitatively  match  those  of  Herrmann  and  Cox  and  that  the  numerical  algorithm  used  to 
implement  the  bond  model  perform  adequately. 
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Figure  19.  -  Shear  Lag  Problem  and  Finite  Element  Model 
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Figure  20.  -  Shear  Lag  Problem  Results 


8.2.1  Simulation  of  Malvar's  Test.  In  Malvar's  (1991)  concrete  bond  tests,  confinement 
pressure  was  applied  to  a  cylinder  of  concrete  surroimding  a  single  reinforcing  bar  which  was 
then  pulled.  The  tests  yielded  bond  shear  and  normal  stresses  as  a  function  of  slip.  This  test  was 
modeled  using  the  driver  program  and  the  simplified  bond  model  described  in  Sections  6  and  7. 
Calculated  bond  shear  stress  versus  slip  and  dilation  are  shown  in  Figures  21  and  22.  These  can 
be  compared  to  calculations  in  Figures  32  and  37  of  Herrmann  and  Cox  (1992).  A  comparison  of 
the  plots  indicates  that  the  simplified  bond  model  and  material  algorithm  produce  results  that  are 
qualitatively  similar  to  those  of  Herrmann  and  Cox. 

8.2.2  Simulation  of  Gambarova's  Test.  The  Gambarova  bond  test  was  somewhat 
similar  to  that  of  Malvar  and  was  also  used  by  Herrmann  and  Cox  in  the  formulation  of  their 
bond  model.  In  this  test,  a  nonsymmetric  reinforced  prism  of  concrete  was  precracked  or  presplit 
parallel  to  the  reinforcement.  The  reinforcement  was  pulled  while  normal  forces,  confinement 
stresses,  were  applied  to  maintain  the  selected  crack  width.  For  more  details  of  the  test  setup 
refer  to  Section  5.3  of  Herrmann  and  Cox  or  to  Gambarova,  Rosati,  and  Zasso  (1989). 

This  test  was  simulated  using  the  driver  program  and  the  simplified  bond  model.  Results 
are  presented  in  the  form  of  shear  and  normal  stress  plots  as  functions  of  slip  in  Figures  23  and 
24,  which  can  be  compared  to  Herrmann  and  Cox's  Figures  29  and  30.  Once  again  these  results 
are  deemed  acceptable  as  the  simplified  bond  model  reproduces  the  important  characteristics  of 
the  response. 

8.3  Structural  Tests 

Having  verified  that  the  plasticity  bond  model  and  algorithm  performed  adequately,  they 
were  implemented  in  the  nonlinear  composite  finite  element  code  as  described  in  Section  7. 
Results  of  a  series  of  tests  on  this  element  for  a  variety  of  reinforced-concrete  structures  follow. 

8.3.1  Simulation  of  Bond  Tension  Tests.  A  variety  of  researchers  have  tested 
reinforced  prisms  of  concrete  to  aid  in  understanding  bond  effects  on  a  reinforced-concrete 
structure.  Houde  and  Mirza  (1974)  have  performed  tests  that  attempt  to  simulate  the  conditions 
in  the  constant  moment  region  of  a  reinforced  beam  between  tensile  cracks.  Houde  and  Mirza 
instrumented  the  reinforcement,  thus  generating  good  data  with  which  to  test  the  bond  model  and 
composite  element  implementation.  The  test  specimen  is  shown  in  Figure  25.  The  single  bar 
reinforcement  was  instrumented  with  strain  gages  in  a  cavity  machined  internally  into  the  bar. 
Loads  were  applied  to  the  bar  in  a  tensile  test  machine.  The  strain  gages  allowed  the  calculation 
of  the  bar  force  distribution  as  the  load  was  transferred  into  the  concrete  prism  through  the  bond. 
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Figure  22.  -  Malvar  Test  Simulatioi^Shear  Stress  vs.  Dilation 
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Figure  23. 
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Figure  24.  -  Gambarova  Test  Simulation-Normal  Stress  vs.  Slip  Strain 
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Figure  25,  -  Houde  and  Mirza  Bond  Tension  Test 


As  previously  noted,  the  composite  finite  element  is  restricted  to  plane  strain  behavior 
which  is  better  suited  to  a  reinforced  slab  than  a  prism.  The  results,  however,  can  still  give  a 
qualitative  feel  for  the  bond  model  and  finite  element  implementation  performance.  The 
composite  finite  element  model  and  boundary  conditions  for  this  test  are  shown  in  Figure  26.  To 
represent  the  hollowed  bar,  the  fiber  radius  and  volume  fraction  were  based  on  the  actual  bar 
diameter  in  order  to  correctly  model  the  interface  in  the  embedded  micro-mechanics  model.  The 
lower  bar  stiffness  resulting  from  the  cavity  was  accounted  for  by  using  an  artificially-reduced 
steel  modulus.  It  is  important  to  note  the  simplicity  of  the  finite  element  model  for  this  test  case. 
Displacements  were  prescribed  to  the  fiber  degree  of  freedom  simulating  the  displacement 
controlled  test  performed  by  Houde  and  Mirza.  The  model  shown  is  the  last  in  a  series  of  three 
models  in  which  the  mesh  size  was  decreased  to  check  for  convergence.  Convergence  was 
obtained  with  one-half  of  the  elements  shown. 

Resulting  predicted  bar  strains,  as  a  function  of  distance  along  the  specimen,  are  plotted 
for  a  variety  of  bar  load  levels  in  Figure  27.  Similar  Houde  and  Mirza  experimental  data  are 
given  in  Figure  28.  The  model  is  able  to  reproduce  the  qualitative  nature  of  the  experimental 
results.  It  is  not  expected  to  match  the  quantitative  results  since  the  model  is  plane  strain  as 
mentioned  earlier.  It  should  also  be  noted  that  bar  lug  spacing  and  angle  which  are  input 
parameters  for  the  bond  model  were  not  given  by  Houde  and  Mirza  and  are  assumed  to  be  similar 
to  those  given  by  other  investigators  who  used  bars  of  similar  diameters. 

The  results  are  very  good  in  light  of  the  known  limitations.  The  shape  of  the  curves  are 
similar  and  the  strain  magnitudes  are  quite  close  for  common  load  levels.  Note  that  the 
composite  element  predicts  detailed  bond  stresses  with  a  very  simple  finite  element  mesh  and 
only  basic  material  parameters. 

8.3.2  Simulation  of  Krefeld  and  Thurston  Cracked-Reinforced  Beam.  The  next 
reinforced  structure  analyzed  was  a  cracked-reinforced  beam.  Krefeld  and  Thurston  (1966) 
tested  a  variety  of  precracked  beams  to  determine  how  the  reinforcing  steel  contributes  to  shear 
resistance.  The  composite  element  as  implemented  in  this  work  cannot  model  the  reinforcement 
dowel  action  that  Krefeld  and  Thurston  were  investigating.  To  model  this  effect,  it  would  appear 
to  be  necessary  to  account  for  finite  deformations  and  the  resulting  finite  length  of  reinforcement 
that  pulls  from  the  surrounding  concrete.  This  essentially  forms  a  beam  spanning  the  opened 
crack.  In  addition,  the  bending  and  shear  of  the  reinforcement  would  need  to  be  accounted  for. 
Krefeld  and  Thurston's  experimental  results  can  still  be  used  to  assess  the  performance  of  the 
composite  element. 
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Figure  26.  -  Houde  and  Mirza  Bond  Tension  Test  Finite  Element  Model 
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Krefeld  and  Thurston  ran  a  series  of  cracked  beam  tests.  Their  test,  titled  DAIO,  was 
selected  for  modeling  because  the  instrumentation  included  strain  gages  located  on  the 
reinforcement  in  the  crack.  The  test  setup  is  shown  in  Figure  29.  The  corresponding  finite 
element  mesh  and  boundary  conditions  are  shown  in  Figure  30. 

The  model  employs  three  different  elements.  The  bottom  three  rows  are  composite 
elements  while  the  upper  elements  are  normal  isotropic  material  elements.  The  quadrilateral 
elements  in  the  crack  are  really  truss  elements  which  have  stiffness  limited  to  the  axial  bar 
response.  No  shear  or  bending  stiffness  was  included.  As  mentioned  earlier,  the  present 
development  is  restricted  to  two-dimensional  plane  strain  and  would  best  model  a  cracked- 
reinforced  slab.  Results  again  are  judged  qualitatively  rather  than  quantitatively  because  of  these 
limitations.  (Note  that  a  plane  stress  implementation  could  be  developed  using  the  three- 
dimensional  material  relations  given  in  Equation  4.1 1 1.  Also  a  three-dimensional  brick  element 
could  be  developed  by  adding  slip  and  dilation  degrees  of  freedom  to  the  three  global 
displacement  degrees  of  freedom  of  a  standard  8-node  hex  element.) 

The  plane  strain  model  was  loaded  by  prescribing  a  vertical  displacement  at  the  beam's 
center.  The  resulting  load  at  the  point  was  accumulated  from  the  global  residual  before 
implementation  of  the  essential  boundary  conditions.  Results  in  the  form  of  a  deformed  mesh 
plot  are  shown  in  Figure  31.  Results  are  consistent  with  the  boundary  conditions  and  the  loss  of 
structure  stiffness  as  a  consequence  of  the  crack.  Figure  32  is  a  collection  of  gray-scale  stress 
contour  plots  for  this  model.  Figure  32a  is  the  matrix  or  concrete  axial  stress.  The  concrete 
stress  goes  to  zero  in  the  vicinity  of  the  crack  as  expected.  Figure  32d  depicts  the  fiber,  or  bar 
axial  stress,  while  Figures  32e  and  32f  show  the  bond  stresses  for  the  composite  elements.  It 
should  be  noted  that  the  contour  plots  were  done  using  averaging  across-neighboring  elements, 
which  affects  the  accuracy  of  stresses  at  the  interface  between  two  material  types.  A  study  of 
these  results  confirms  the  expected  behavior  of  a  beam  cracked  in  this  fashion. 

Figure  33  gives  experimental  results  from  Krefeld  and  Thurston  (1966)  in  which  the  force 
in  the  reinforcing  bars  is  plotted  versus  the  applied  machine  load.  Figure  34  is  a  similar  plot 
using  the  composite  finite  element  results.  The  model  is  not  able  to  predict  the  nonlinearities 
present  at  the  higher  load  levels.  These  nonlinearities  are  assumed  to  be  a  result  of  concrete 
inelastic  response  due  to  crushing  and  cracking  and  bar  yield.  Recall  that  the  composite 
element’s  nonlinear  response  is  limited  to  that  of  the  bond  and,  thus,  carmot  be  expected  to 
predict  such  behavior.  The  good  quantitative  agreement  at  the  lower  loads  again  verifies  the 
composite  element  approach.  Figure  35  presents  more  finite  element  results  where  the  bar  and 
machine  loads  are  plotted  as  a  function  of  the  mid-span  displacement.  Note  that,  in  the  case  of 
the  machine  load,  one  can  see  the  effects  of  the  bond  nonlinearity  as  the  structure  softens  as  a 
result  of  slip. 
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Figure  29.  -  Krefeld  and  Thurston  Dowel  Test  DA  10 
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Figure  30.  -  Krefeld  and  Thurston  Beam  -  Finite  Element  Plane  Strain  Model 


Figure  31.  -  Krefeld  and  Thurston  Model,  Deformed  Mesh  Results 
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Figure  32a.  -  Concrete  Matrix  x-Direction  Stress 
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Figure  32b.  -  R/C  Composite  y-Direction  Stress 
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Figure  32.  -  Krefeld  and  Thurston  Beam  Simulation 
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Figure  32c.  -  R/C  Composite  xy-Shear  Stress 
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Figure  32d.  -  Bar  x-Direction  Stress 
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Figure  32e.  -  Bond  Shear  Stress 
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Figure  32f.  -  Bond  Normal  Stress 
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Figure  34.  -  Krefeld  and  Thurston  Simulation;  Bar  vs.  Machine  Load  Finite  Element  Results 
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This  model  also  served  as  a  rigorous  test  of  the  bond  material  algorithm.  The  initial 
results  of  this  analysis  motivated  important  needed  changes  and  yielded  lessons  on 
implementation  and  usage.  It  was  stated  previously  that  the  study  was  limited  to  monotonic 
loading.  In  initial  attempts  to  analyze  this  test,  it  was  observed  that  monotonic  structural  loading 
can  still  result  in  cyclic  loading  at  the  element  level.  As  slip  evolves,  a  region  may  go  from  slip 
in  one  direction  to  that  of  the  other  as  the  concrete  flexural  crack  edge  effects  reach  further  into 
the  structure.  The  bond  model  was  modified  to  handle  this  simple  one-cycle  behavior.  A  very 
simple  implementation,  which  assumed  that  slip  in  either  direction  results  in  damage  and 
evolution  of  the  yield  surface,  was  conceived.  While  this  is  known  to  be  incorrect  (see  the  cyclic 
model  of  Herrmann  and  Cox,  1992),  it  was  deemed  adequate  for  this  case  where  the  amount  of 
slip  in  the  initial  direction  is  very  small  and  does  not  contribute  significantly  to  the  damage. 

Another  important  observation  made  in  the  initial  attempts  to  analyze  this  test  was  that 
the  evaluation  of  the  bond  material  algorithm  for  some  elements  is  difficult  for  a  real  structure 
where  there  are  regions  of  low  confinement  or  near  zero  bond  normal  stress.  The  load  had  to  be 
applied  very  slowly  in  the  incremental  scheme  for  the  entire  loading  sequence  to  prevent 
numerical  instabilities.  However,  there  are  also  limits  as  to  how  small  the  load  steps  could  be 
made  because  of  numerical  instability,  possibly  as  a  result  of  loss  of  precision.  The  distinction 
between  too  big  or  too  small  a  load  step  was  subtle  in  some  cases,  making  it  difficult  to  know 
how  much  load  or  displacement  to  apply  in  a  given  increment.  The  complete  resolution  of  this 
problem  will  require  additional  study  in  future  work. 

8.3.3  Simulation  of  Bressler  and  Scordelis  Reinforced-Concrete  Beam.  A  second 
reinforced  beam  structure  was  analyzed.  A  series  of  beam  tests  were  performed  by  Bressler  and 
Scordelis  (1963)  to  study  the  shear  strength  of  reinforced-concrete  beams.  The  series  OA-1  was 
selected  to  model.  The  beam  structure  and  test  setup  are  shown  in  Figure  36.  The  beam  was 
initially  uncracked  but  the  experimenters  tracked  and  recorded  flexural  and  diagonal  tension 
cracks  as  they  incrementally  loaded  the  beam.  Once  again,  the  present  composite  element  cannot 
be  expected  to  model  this  complicated  evolution  and  behavior.  To  do  this  would  require  as  a 
minimum  an  inelastic  concrete  material  model  coupled  with  the  ability  to  spawn  concrete  cracks 
in  both  composite  and  concrete  elements. 

However,  it  is  still  possible  to  test  the  present  element’s  performance  by  assuming  initial 
flexviral  cracks  derived  from  Bressler  and  Scordelis’  reported  crack  patterns.  This  model  is 
shown  in  Figure  37.  The  model  is  plane  strain  which  again  is  a  limitation.  It  employs  two  rows 
of  composite  elements  to  represent  the  bottom  portion  of  the  beam  where  the  four  reinforcing 
bars  were  located.  There  are  thin,  bar-only  quad  elements  between  the  composite  elements  at  the 
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assumed  crack  locations.  The  1 3  cracks,  for  the  full  beam,  extend  through  the  bottom  three  rows 
of  elements.  Additionally,  preliminary  one-crack  and  seven-crack  models  were  also  analyzed. 


Model  Data: 
f 'c  =  22.6  N/mm2 

p  =  .044  for  bottom  7.5  in  of  beam 
Econcrete=37.5  N/mm2 


Figure  36.  -  Dressier  and  Scordelis  Beam  Test  OA-1 
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Figure  37.  -  Bressler  and  Scordelis  Beam  -  Finite  Element  Model 


The  model  was  analyzed.  The  resulting  deformed  mesh  plot  is  shown  in  Figure  38.  The 
matrix  cracks  open  as  slip  occurs.  Stress  contour  plots  are  given  in  Figure  39.  Again  axial 
tensile  stress  cannot  develop  in  the  concrete  because  of  the  cracks  and,  thus,  flows  around  them 
as  can  be  seen  in  Figure  39a.  Figure  39c  shows  the  reversal  in  shear  stress  across  the  crack  in  the 
reinforced  concrete.  Figure  39d  shows  the  gradient  in  bar  tensile  stress  due  to  bond  along  the  bar 
length.  Figure  39e  emphasizes  the  bond  shear  stress  behavior  between  cracks.  The  contours 
show  perfectly  the  locations  of  the  cracks.  Across  a  crack,  the  sign  of  the  shear  stress  switches  as 
it  should.  Approximately  halfway  between  cracks  the  bond  shear  stress  goes  to  zero.  This  can 
also  be  observed  by  reviewing  the  bond  normal  stress  plotted  in  Figure  39f.  In  this  case,  the 
confinement  pressure  peaks  at  the  crack  and  drops  off  to  zero.  The  composite  element  was  able 
to  capture  this  detailed  bond  behavior  provided  the  specified  load  or  displacement  was  applied 
slowly  enough. 

Figure  40  compares  experimental  and  analytical  load-deflection  behavior.  Once  again, 
one  cannot  expect  close  quantitative  agreement  in  light  of  all  the  restrictions  mentioned 
previously  and  the  fact  that  the  model  is  plane  strain  and  starts  with  13  cracks.  The  results  are 
judged  to  be  quite  good  considering  all  the  limitations  and  approximations. 


9.  CONCLUSIONS 

The  results  of  this  research  clearly  demonstrate  the  feasibility  of  incorporating  a  bond 
model,  that  includes  both  bond  slippage  and  dilation,  into  a  composite  representation  and 
accompanying  finite  element  analysis  for  unidirectional  reinforced  materials. 

A  three-dimensional  micro-mechanics  model  of  an  approximate  representative  volume 
for  a  imidirectional  composite  was  formulated.  The  model  assumes  linear-elastic  properties  for 
the  reinforcing  and  matrix  materials  and  a  plasticity  model  for  the  bond  between  them,  including 
bond  dilation.  The  resulting  composite  model  was  successfully  incorporated  into  a  two- 
dimensional  composite  finite  element  analysis. 

Several  sample  analyses  were  successfully  performed  for  simple  reinforced-concrete 
tests.  The  analysis  of  the  bond  tension  test  served  to  demonstrate  best  the  current  capability  of 
the  composite  element.  The  mesh  was  a  simple  row  of  quad  elements.  The  results  that  were 
obtained  were  quite  detailed.  For  example,  the  fiber  or  bar  stress  was  determined  and  the  bond 
stresses  were  predicted.  Such  results  are  normally  only  associated  with  a  detailed,  discrete  finite 
element  model  of  the  reinforcement  with  an  interface  element  and  provision  for  bond-slip 
behavior.  The  composite  element  approach  has  been  shown  to  be  a  powerful,  efficient  and 
capable  method  of  modeling  bond  behavior  in  reinforced  materials. 
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Figure  38.  -  Bressler  and  Scordelis  Beam  -  Finite  Element  Deformed  Mesh 
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Figure  39a.  -  Concrete  Matrix  x-Direction  Stress 
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Figure  39b.  -  R/C  Composite  y-Direction  Stress 
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Figure  39.  -  Dressier  and  Scordelis  Beam  Simulation 
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Figure  39c.  -  R/C  Composite  xy-Shear  Stress 
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Figure  39d.  -  Bar  x-Direction  Stress 
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Figure  39e.  -  Bond  Shear  Stress 
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Figure  39f.  -  Bond  Normal  Stress 
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Figure  40.  -  Load  vs.  Deflection  Results  for  Dressier  and  Scordelis  Beam 


The  element  employs  two  additional  nodal  degrees  of  freedom  to  capture  slippage  of  the 
reinforcement  relative  to  the  matrix  and  the  dilation  of  the  bond  zone.  The  use  of  the  four 
degrees  of  freedom  (two-dimensional)  quad  element  is  simple  compared  to  a  discrete  model 
which  would  require  finely  meshed  patches  of  elements  for  each  constituent  type,  making  it 
impractical  to  model  bond-slip  behavior  at  the  structure  level.  The  composite  element  results  in 
savings  in  both  the  construction  of  the  finite  element  model  and  the  subsequent  bond-slip 
analysis. 

The  virtual  work  micro-mechanics  approach  proved  to  be  an  effective  means  for  deriving 
the  specialized  composite  material  relations  used  in  the  finite  element  scheme.  The  Reduced- 
Newton  material  algorithm,  developed  to  evaluate  the  interface  bond  model,  is  quite  efficient  and 
reasonably  robust. 


10.  RECOMMENDATIONS 

While  much  was  accomplished  in  this  study,  additional  work  is  needed  before  the  finite 
element  analysis  can  be  applied  to  significant  reinforced-concrete  problems.  The  model  must  be 
extended  to  include  inelastic  properties  for  the  concrete  and  possibly  for  the  reinforcing  bars,  and 
some  means  must  be  included  for  modeling  the  concrete  cracking.  For  the  latter  extension,  the 
introduction  of  smeared  cracking  (possibly  by  using  a  generalized  plasticity  model  for  the 
concrete)  would  seem  the  most  appropriate.  The  extension  of  the  finite  element  analysis  to 
generalized  plane  stress  and  three-dimensional  conditions  would  greatly  increase  the  range  of 
reinforced-concrete  problems  that  could  be  addressed. 

A  second  major  thrust  would  involve  the  application  of  the  bond-slip  model  and  finite 
element  analysis  to  fiber-reinforced  composites.  The  key  step  would  be  to  use  whatever 
experimental  data  are  available  to  calibrate  (and  modify  as  necessary)  the  bond  model  for  one  or 
more  composite  systems  of  this  type.  The  finite  element  analysis  could  then  be  used  to  analyze 
any  available  structural  tests  of  fiber-reinforced  components  and  a  judgment  could  be  made  as  to 
the  importance  of  including  bond  slippage  and  dilation  in  the  modeling  of  such  systems. 

Finally,  three  components  of  the  current  composite  model  deserve  further  study;  (1)  the 
inclusion  of  the  a*  variable,  the  average  bond  equilibrium  stress,  in  the  stress  vector  is 
cumbersome  -  initial  attempts  to  find  an  alternative  scheme  for  specifying  interface  equilibrium 
were  not  successful  but  with  additional  investigation  one  might  be  found;  (2)  both  the  theoretical 
and  practical  advantages  and  disadvantages  of  using  element  condensation  to  eliminate  the  bond 
dilation  from  the  global  degrees  of  freedom  need  to  be  addressed  in  future  work;  and  (3)  an 
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automated  global  loading  algorithm  to  eliminate  time-consuming  trial-and-error  selection  of  load 
increment  sizes  is  required  for  robust  applications. 
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